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1.  INTRODUCTION 


1 . 1  Summary 

It  is  known  that  large  space  structures  will  be  subjected  to  thermome¬ 
chanical  loadings  and  environmental  conditions  which  are  likely  to  degrade  the 
constitutive  properties  of  the  structural  materials,  thus  leading  to  possible 
failure  of  these  vehicles.  Therefore,  it  is  desirable  to  develop  new  analytical 
models  which  are  capable  of  accounting  for  these  degraded  properties  so  that  de¬ 
sign  procedures  can  be  improved.  There  are  three  important  aspects  of  such  an 
effort:  selection  and  development  of  constitutive  models  applicable  to  large 
space  structures,  construction  of  analytic  models,  and  experimentation  to  deter¬ 
mine  the  precise  nature  of  the  material  parameters  to  be  utilized  in  the  analyt¬ 
ical  model.  These  three  components  of  the  research  must  be  tied  together  into 
a  single  concise  package  in  order  to  obtain  a  useful  model. 

This  research  project  is  a  three  year  effort  to  develop  an  analytic  model 
capable  of  predicting  the  response  of  space  structures  with  degrading  material 
properties  under  quasi-static  as  well  as  dynamic  cyclic  thermomechanical  loading 
conditions.  This  report  details  the  research  completed  during  the  first  year 
of  AFOSR  contract  no.  F49620-83-C-0067 . 

1 . 2  Statement  of  Work 

A  model  is  being  developed  for  predicting  the  thermomechanical  response  of 
large  space  structures  to  cyclic  transient  temperature  loading  conditions.  The 
research  is  being  conducted  in  the  following  stages: 

1)  selection  and  specialization  of  thermomechanical  constitutive  equations 
to  be  utilized  in  the  analysis  of  large  space  structures; 

2)  construction  (where  necesssary)  of  coupled  energy  balance  equations 
(modified  Fourier  heat  conduction  equations)  applicable  to  the  consti¬ 
tutive  models  selected  in  item  1); 


3)  casting  (where  necessary)  the  resulting  field  laws  into  coupled  and 
uncoupled  variational  principles  suitable  for  use  with  the  finite 
element  method; 

4)  finite  element  discretization  of  the  variational  principles  for 
several  element  types; 

5)  experimentation  to  determine  material  properties  to  be  utilized  in  the 
constitutive  models;  and 

6)  parametric  studies  of  the  quasi-static  and  dynamic  response  of  large 
space  structures  undergoing  thermomechanically  and  environmentally 
degraded  material  properties. 

The  experimental  effort  (discussed  in  5)  is  being  supported  in  part  by  DOD 
equipment  grant  no.  841542  with  equipment  on-line  capability  expected  by  December 
31,  1984.  The  total  research  effort  outlined  above  spans  a  period  of  three  years. 
The  following  section  details  results  obtained  during  the  first  year. 


2.  RESEARCH  COMPLETED  TO  DATE 

2 . 1  Summary  of  Completed  Research 

The  following  tasks  have  been  completed  during  the  first  year  of  research: 

1)  literature  survey; 

2)  selection  of  constitutive  equations  for  thermoviscoplastic  metals  at 
elevated  temperatures  and  polymeric  composites  with  thermomechanical 
load  induced  damage  at  temperatures  below  the  glass  transition  tempera¬ 
ture; 

3)  construction  of  a  coupled  energy  balance  equation  for  thermovisco¬ 
plastic  metals; 

4)  casting  of  field  laws  for  the  materials  discussed  in  2)  into 
one-dimensional  (truss  and  beam)  finite  element  computer  codes  with 
one-  and  two-way  thermomechanical  coupling;  and 

5)  initial  parametric  studies  using  the  models  developed  in  4)  to  determine 
the  thermomechanical  response  of  representative  space  structures  with 
degraded  material  properties. 

The  work  statement  for  the  second  year  is  as  follows: 

1)  refinement  of  constitutive  equations  for  polymeric  composites  based  on 
currently  available  research; 

2)  selection  and  development,  to  include  experimental  characterization, 
of  metal  matrix  composite  constitutive  equations  appropriate  for  large 
space  structures; 

3)  construction  of  coupled  energy  balance  equations  for  metal  matrix 
composites  and  application  of  the  resulting  field  laws  to  truss  finite 
element  codes; 

4)  further  development  of  beam  finite  element  codes  for  viscoplastic  metals, 
polymeric  composites,  and  metal  matrix  composites  to  include  one-  and 
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two-way  thermomechanical  coupling; 

5)  determination  of  the  effect  of  inelastic  heat  generation  on  the  thermo- 
vibratory  response  of  representative  space  structures;  and 

6)  completion  of  further  parametric  studies  using  the  models  developed  in 
step  4) . 

In  addition,  a  DEC  11/23  (LSI  11/23  Plus)  has  been  purchased  with  internal 
university  funding  and  is  now  being  utilized  to  generate  many  of  the  numerical 
results  obtained  on  the  contract. 

Additional  experimental  equipment  is  also  being  purchased  under  DOD  equip¬ 
ment  grant  no.  841542.  This  equipment  will  support  item  2)  of  the  second  year 
work  statement  listed  above.  At  this  time  the  equipment  has  been  ordered  and 
is  expected  to  be  on-line  by  December  31,  1984. 

In  addition  to  the  funding  supplied  by  AFOSR,  the  Texas  Engineering  Experi¬ 
ment  Station,  a  branch  of  the  Texas  A&M  University  System,  has  provided  sub¬ 
stantial  internal  support  for  the  research  detailed  herein.  This  support  has 
been  utilized  to  support  additional  graduate  students  and  for  travel.  In  addi¬ 
tion,  Dr.  Michael  S.  Pilant  of  the  Mathematics  Department  at  TAMU  is  being 
supported  for  one-half  month  per  year  to  study  computational  algorithms  for 
solving  thermoviscoplastic  field  problems. 

Finally,  it  should  be  pointed  out  that  a  material  property  degradation 
experiment  is  currently  underway  on  the  NASA  LDEF  under  the  direction  of  Dr. 
Allen.  Effects  of  long  term  space  exposure  on  material  properties  of  poly¬ 
ethylene  will  become  available  during  the  second  contract  year. 

The  following  sections  detail  the  results  obtained  during  the  first  con¬ 
tractual  year. 
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2.2  Literature  Surve 


A  detailed  literature  survey  has  been  completed  as  part  of  the  first  year 
effort.  This  report,  entitled  "Large  Space  Structures  Technology:  A  Literature 
Survey",  is  included  as  Appendix  6.3.  Briefly,  the  report  details  recent  ad¬ 
vances  in  the  areas  of  materials,  structural  solution  techniques,  damping,  and 
preliminary  design  and  experiment.  The  results  of  this  survey  indicate  that 
very  little  research  is  available  on  the  effect  of  material  property  degrada¬ 
tion  on  large  space  structural  response. 


2.3  Selection  of  Constitutive  Equations 

Candidate  material  models  have  been  selected  for  metals  at  elevated  tem¬ 
peratures  and  polymeric  composites  below  the  glass  transition  temperature. 

These  are  detailed  below. 

2.3.1  Metals  at  elevated  temperatures  are  currently  modeled  using  continuum 
mechanics  with  internal  state  variables  (ISV's)  [l-5],  wherein  the  stress-strain 
relation  is  of  the  form  (for  infinitesimal  strains): 


I  T 

a .  .  -  D  (e.  -  e.  -  e.  ) 

ij  ijkj,  k.2,  kd  k£ 


(1) 


where  a.,  is  the  stress  tensor,  e,  is  the  strain  tensor,  D...  is  the  linear 
ij  kd  xjkd 

I  T 

elastic  modulus  tensor,  e^£  is  the  inelastic  strain  tensor,  and  e  is  the 

thermal  strain  tensor.  In  addition. 


ij 


fi(Ew  T’  V 


(2) 


and 


W  ^ 


ij 


f  ^  ( £  T  or  } 
rij^  kd* 


(3) 


where  f  and  f.  are  appropriate  functions  of  state,  T  is  the  temperature,  and 


k.. 


are  a  set  of  second  order  tensor  valued  internal  state  variables  modeling 


dislocation  arrangement,  dislocation  density,  intergranular  damage,  etc. 
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Although  it  has  been  demonstrated  that  numerous  models  fall  within  the 
above  framework  [5],  the  special  cases  of  equations  (1)  through  (3)  utilized 
thus  far  are  a  classical  plasticity  model  developed  by  Allen  and  Haisler  [6,  7] 
(see  Appendix  6.2),  a  single  internal  state  variable  viscoplastic  model  developed 
by  Cernocky  and  Krempl  [8,  9]  (see  Appendix  6.1),  and  a  two  internal  state  vari¬ 
able  viscoplastic  model  developed  by  Bodner,  et  al.  [10-13]  (see  Appendices  6.1 
and  6.5).  It  is  emphasized,  however,  that  the  algorithms  developed  under  this 
contract  can  be  utilized  with  any  model  capable  of  formulation  according  to 
equations  (1)  through  (3). 

2.3.2  Polymeric  composites  are  modeled  using  internal  state  variable  theory 
as  well.  However,  in  this  case  the  ISV's  are  assumed  to  be  vector-valued  and  to 
enter  only  through  the  modulus  tensor,  that  is: 


a .  . 

iJ 


=  c:..  ( e . .  —  e :.) 

ijk£  ij  ij' 


(A) 


where  C'.  , ,  the  effective  modulus  tensor,  is  given  by 
1J  Kx 


Cl  =  C.  -  PPq  .  cxPaq 
13  13  kx  mnkicij  m  n 


where  aP  are  a  set  of  vector-valued  internal  state  variables  [lA]  describing 
transverse  cracking,  interlaminar  delamination,  fiber  breakage,  etc.,  and 
modeled  by 

“i  ■  ai<Ekt-  t-  •  <6> 


At  low  homologous  temperatures  these  materials  are  assumed  to  rate  insensi¬ 
tive  so  that  the  above  model  will  result  in  quasi-elastic  equations  in  which 
inelasticity  is  reflected  only  through  the  slowly  degrading  modulus  tensor. 
Experimental  evidence  [15,  16]  indicates  that  the  time  scale  for  degradation 


of  C 


ijk£ 


is  very  long  compared  to  the  frequencies  and  mode  shapes  of  representat ive 


structures.  It  is  therefore  sufficient  for  many  space  structural  applications 
to  treat  equations  (4)  in  the  degraded  state  only.  This  approach  is  further 
detailed  in  Appendix  6.4. 

2.3.3  Metal  matrix  composites  are  currently  being  purchased  from  ARCO 
Metals,  Greenville,  S.D.,  and  models  for  these  materials  will  be  developed 
during  the  second  and  third  years  of  the  contract. 


2.4  Coupled  Energy  Balance  Laws 

The  energy  balance  law  for  thermomechanically  coupled  media  of  the  type 
described  in  section  2.3.1  has  been  constructed  [ 1 7]  (see  also  Appendix  6.1). 
This  equation  can  be  utilized  to  predict  temperature  rise  in  a  thermovisco¬ 
plastic  medium  subjected  to  cyclic  mechanical  loading.  This  equation  is  in 
general  a  statement  of  conservation  of  energy  and  represents  a  modification  of 


the  Fourier  heat  conduction  equation  given  by 


D..,  (e  -a  +a,  T  )a1..+D...  a. .a,  TT  -  D  a  .  .  Te  -pC  T  +  q.,.=pr 
ljkS.  kl  lk Si  kd  R  lij  xjki  lj  k£  ijk2  i]  ki  v  j  j 


where  a  is  the  internal  state  variable  modeling  drag  stress,  a  is  the 

J.  iv  jG  iC  X. 

thermal  conductivity  tensor,  p  is  the  mass  density,  Cv  is  the  specific  heat, 
q^  is  the  heat  flux  vector,  and  r  is  the  specific  heat  supply. 

The  above  result  has  been  utilized  to  predict  the  thermomechanical  response 
of  a  single  insulated  truss  element  to  cyclic  mechanical  loading  (see  Appendix 
6.1).  As  shown  in  Figs.  1  and  2,  substantial  temperature  rise  (approximately 
3.7°C)  is  predicted  for  each  cycle. 

It  is  now  thought  that  in  polymeric  composites  about  95%  of  the  strain  energy 
lost  to  inelastic  deformations  is  expended  in  the  creation  of  internal  surfaces 
called  damage.  It  is  therefore  not  necessary  to  construct  two-way  coupled  energy 
balance  laws  for  these  materials  and  the  classical  Fourier  heat  conduction  equa¬ 
tion  is  adequate  for  modeling  the  temperature  field.  Therefore,  the  models  de¬ 
veloped  herein  utilize  only  one-way  coupling  for  polymeric  composite  media:  that 
is,  the  temperature  field  affects  the  displacement  field  but  not  vice  versa. 
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2.5  Finite  Element  Computer  Codes 

Five  finite  element  computer  codes  have  been  developed  and/or  refined  dur¬ 
ing  the  first  year  of  research.  These  are: 

1)  one-dimensional  truss  type  heat  conduction  code  for  predicting 
time-dependent  axial  temperature  variations  in  truss  structures; 

2)  one-dimensional  truss  code  for  predicting  dynamic  response  of  truss 
structures  with  quasi-elastic  materials  properties  to  thermomechanical 
inputs  (see  Appendix  6.4); 

3)  one-dimensional  beam  code  for  predicting  dynamic  response  of  beam 
structures  with  quasi-elastic  material  properties  to  thermomechanical 
inputs  (see  Appendix  6.4); 

4)  one-dimensional  fully  two-way  coupled  truss  code  for  predicting  the 
quasi-static  response  of  truss  structures  with  thermoviscoplastic 
material  properties  to  thermomechanical  inputs  (see  Appendix  6.5); 
and 

5)  two-dimensional  continuum  code  for  predicting  quasi-static  response  of 
elastic-plastic  media  to  mechanical  inputs  (see  Appendix  6.2). 

Algorithms  for  constructing  elastic  versions  of  items  1)  through  3)  are 
in  the  research  literature  and  will  not  be  discussed  further  herein.  Details 
of  the  quasi-elastic  material  properties  are  discussed  in  Appendix  6.4.  Results 
of  the  analytical  models  are  significant  in  that  they  show  substantial  changes 
can  occur  in  the  dynamic  response  of  representative  space  structures  due  to 
small  changes  in  the  material  properties  of  the  structural  elements.  A  brief 
review  of  these  results  will  be  given  in  the  following  section. 

The  fully  two-way  coupled  truss  code  has  been  constructed  at  this  time,  as 
detailed  in  Appendix  6.5.  However,  numerical  results  have  not  yet  been  obtained. 
These  should  be  forthcoming  in  the  near  future. 
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Finally,  the  numerical  efficiency  of  a  previously  developed  elasto- 
plasticity  code  has  been  significantly  increased,  as  detailed  in  Appendix  6.2. 


The  codes  listed  above  will  be  utilized  and  further  refined  by  the  inclu¬ 
sion  of  more  advanced  material  models  during  the  second  and  third  years  of  the 
research  effort. 

2 . 6  Initial  Parametric  Studies 

For  purposes  of  illustrating  the  capability  of  the  currently  developed 
models,  several  sample  problems  have  been  provided  for  space  structures  with 
realistic  degradation  of  material  properties  (see  Appendix  6.4).  As  an  example, 
consider  the  truss  structure  shown  in  Fig.  3.  This  beam-like  structure  is  canti¬ 
levered  at  one  end  to  simulate  an  antenna  boom.  The  structure  is  sixty  feet 
long  with  bays  ten  feet  long  by  three  feet  wide.  The  structure  is  constructed 
from  graphite-epoxy  composite  material  with  a  quasi-isotropic  ply  layup.  Experi¬ 
mental  research  indicates  that  the  material  may  undergo  up  to  15  percent  loss  in 
stiffness  due  to  cyclic  thermomechanical  fatigue  which  causes  a  variety  of  damage 
modes  in  the  structure.  Additional  loss  of  stiffness  may  be  attributed  to  ele¬ 
vated  temperature  and  chemical  changes  due  to  solar  radiation  and  other  environ¬ 
mental  effects.  In  this  model  the  properties  are  degraded  spacially  on  an 
element  by  element  basis  as  a  function  of  the  stress  history  in  the  structure 
induced  by  long  term  thermomechanical  cyclic  loading.  Stress  amplitudes  were 
obtained  by  using  displacements  corresponding  to  the  first  modal  shape  and  the 
degraded  properties  were  computed  by  assuming  a  linear  damage  law  based  on  peak 
stress  amplitude.  Because  the  boom  is  fixed  on  one  end,  the  stresses  are  high- 


this  spacially  induced  damage.  Figure  5  shows  that  the  shape  of  the  first  mode 
undergoes  no  appreciable  change  in  shape  as  the  damage  occurs,  which  is  due  to 
the  fact  that  the  first  mode  is  a  symmetric  mode.  The  second  mode  shape,  how¬ 
ever,  shown  in  Fig.  6,  undergoes  a  substantial  change  in  shape.  This  result 
indicates  that  active  control  mechanisms  which  may  be  placed  according  to  the 
original  undamaged  mode  shapes  may  not  be  capable  of  controlling  all  modes  as 
the  structure  dynamic  response  changes  over  long  time  periods.  The  control 
designer  must  be  cognizant  of  these  modal  changes  if  he  is  to  design  a  workable 
control  system.  Adaptive  or  robust  control  systems  will  be  required. 

Numerous  other  results  are  presented  in  Appendix  6.4.  For  example,  it  is 
shown  that  if  the  material  properties  degrade  on  only  onde  side  of  the  struc¬ 
ture,  as  might  happen  due  to  solar  radiation,  the  mode  shapes  are  completely 
changed.  Additional  results  are  presented  for  the  transient  responses  of  truss 
and  frame  structures  subjected  to  thruster  loads  and  to  cyclic  solar  radiation 
(heating)  loads. 

Further  refinement  of  the  material  models  is  expected  to  improve  the  capa¬ 
bility  of  the  analytic  model  to  predict  the  effects  of  material  inelasticity  on 
dynamic  response  of  large  space  structures. 

In  conclusion,  these  preliminary  results  seem  to  indicate  that  small  changes 
(or  errors)  in  material  properties  as  they  change  or  degrade  due  to  fatigue, 
damage,  etc.  may  produce  significant  changes  in  predicted  dynamic,  frequency 
and  modal  response.  Correspondingly,  this  affects  our  ability  to  design  effec¬ 
tive  control  systems  and  places  an  even  greater  burden  on  the  control  designer 
to  develop  systems  which  account  for  these  structural  changes.  It  is  clear  that 
an  understanding  of  material  behavior  in  space  environments  and  its  impact  on 
structural  response  is  very  important  to  successful  design  and  development  of 
large  space  structures. 
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A  Prediction  of  Heat  Generation  in  a 
Thennoviscoplastic  Uniaxial  Bar 

by 

David  H.  Allen 

Abstract 

A  thermodynamic  model  is  presented  for  predicting  the  thermomechanical 
response,  including  temperature  change,  in  a  uniaxial  bar  composed  of  a 
thermoviscoplastic  metallic  medium.  The  model  is  constructed  using  the 
concept  of  internal  state  variables,  and  it  is  shown  that  this  general  frame¬ 
work  is  capable  of  encompassing  several  constitutive  models  currently  used 
to  predict  the  response  of  rate  sensitive  metals  in  the  inelastic  range. 

Results  are  obtained  for  monotonic  loading  which  agree  with  predicted  re¬ 
sults  previously  obtained  by  Cernocky  and  Krempl  for  mild  steel  at  room 
temperature.  The  model  is  then  utilized  in  conjunction  with  Bodner  and 
Partora's  constitutive  equations  to  predict  temperature  change  in  Inconel 
(IN)  100  subjected  to  both  monotonic  and  cyclic  loading  at  1005 °K(1350°F) . 

Introduction 

It  has  long  been  known  that  mechanical  and  thermodynamic  coupling  ex¬ 
ists  in  solid  bodies  [1,2].  However,  in  elastic  bodies  this  coupling  is 
negligible  except  when  mass  inertia  is  not  negligible  due  to  flux  of  heat  generated 
through  the  boundary  of  the  body  [3].  However,  in  thermoviscoplastic  metals 
the  conversion  of  mechanical  energy  to  heat  may  be  singificant  even  under 
non-inertial  conditions,  especially  since  material  properties  become  ex¬ 
tremely  temperature  sensitive  in  the  inelastic  range  of  response  [4-11]. 

Similar  research  has  been  performed  on  non-metallic  media  [12-15]. 


General  concinuum  mechanics  models  have  been  formulated  for  broad  classes 
of  materials  [16-19].  However,  to  this  author's  knowledge  only  recently 
has  attention  been  paid  to  the  coupled  heat  conduction  equation  for  thermo¬ 
viscoplastic  metals  [11,20].  Recently,  Cernocky  and  Krempl  [11]  have  pro¬ 
posed  a  model  for  predicting  the  temperature  rise  in  a  class  of  thermovis¬ 
coplastic  metals,  with  special  emphasis  on  test  coupons  subjected  to  either 
homogeneous  uniaxial  or  torsion  loadings.  In  this  paper  an  alternative  ap¬ 
proach  to  that  proposed  in  [11]  is  discussed.  This  method  uses  the  thermo¬ 
dynamics  with  internal  state  variables  originally  reported  in  [17]  and  dis¬ 
cussed  elsewhere  in  detail  for  metals  [18,21,22],  with  development  of  the 
multidimensional  coupled  heat-conduction  equation  in  [20]. 

The  research  herein  is  presented  in  three  parts:  field  formulation 
in  one-dimensional  form;  development  of  the  governing  equations  from  addi¬ 
tional  constitutive  assumptions;  and  numerical  results  for  selected  problems. 

Thermodynamics  of  a  Uniaxial  Thermoviscoplastic  Bar 

Consider  a  slender  bar  which  is  subjected  to  a  homogeneously  applied 
deformation  field  such  that  the  resulting  stress  field  is  everywhere  uniax¬ 
ial  in  the  x^=x  coordinate  direction,  as  shown  in  Fig.  1.  Rigor  would  re¬ 
quire  that  the  possibility  of  finite  deformations  be  considered.  However 
this  condition  is  covered  in  detail  elsewhere  [17,18,20,21,22],  and  for 


purposes  of  simplicity  only  infinitesimal  deformations  will  be  considered 


herein. 

For  notational  simplicity, 

then,  the  observable  mechanical 

state 

lx 

variables 

are 

u  = 

u^  *  deformation  field, 

(l) 

£ 

= 

*  infinitesimal  strain 

field,  and 

(2) 

tv 

,  = 

7.  ,  ■  stress  f ield 

L  1 

(3) 

Cv 

undeformed 


deformed 


Although  transverse  components  of  deformation  and  strain  may  occur,  it  is 


assumed  that  they  are  not  necessary  to  characterize  the  uniaxial  stress  j. 

The  mechanical  state  variables  (1)  through  (3)  are  adjoined  with  the 
thermodynamic  state  variables: 

e  =  internal  energy  per  unit  mass;  (4) 

r  =  heat  supply  per  unit  mass;  (5) 

s  =  entropy  per  unit  mass,  (6) 

T  =  absolute  temperature;  and  (7) 

q  =  =  heat  flux  in  the  x^  coordinate  direction  ,  (8) 

where  it  is  assumed  in  (8)  that  the  bar  is  isotropic  and  long  and  slender 

with  perfectly  longitudinal  boundaries  so  that  the  heat  flux  is  one-dimen¬ 
sional. 


In  accordance  with  the  theory  of  internal  state  variables  [17],  ob¬ 
servable  state  variables  (1)  through  (8)  are  now  supplemented  with  inter¬ 
nal  state  variable  growth  laws  in  order  to  characterize  the  state  of  inel¬ 
astic  bodies: 

=  scalar  valued  internal  state  variables,  k  ■  i,  ...  ,  n;  (9) 

where  n  is  the  number  of  internal  state  variables  required  to  characterize 
the  state  of  the  body.  The  precise  nature  of  (9)  will  be  discussed,  later. 

Parameters  (1)  through  (9)  are  assumed  to  be  functions  of  space  (x) 
and  time  (t) ,  and  are  assumed  to  be  sufficient  to  describe  the  uniaxial 
state  of  the  bar  at  all  times.  These  parameters  are  constrained  by: 

a)  the  conservation  of  momentum, 

if  -  °  .  ao) 

where  inertial  effects  and  the  body  force  are  assumed  to  be  negligible; 

b)  the  strain-displacement  relation, 


3  U 
?:< 


(11) 


c)  Che  balance  of  energy, 


(12) 


(13) 


where  '(  is  called  che  internal  entropy  race  per  unit  mass. 

As  detailed  by  Coleman  &  Gurtin  [17],  equations  (10)  through  (13)  are 
now  supplemented  with  the  following  constitutive  assumptions: 


7  -  He,  T,  3T/3x,  c^)  ,  (14) 

e  *  e(7,  T,  3T/3x,  ct^)  ;  (15) 

s  =  s  ( € ,  T,  3T/3x,  c^)  ;  (16) 

q  -  q (c ,  T,  3T73x,  o^)  ;  and  (17) 

*  \(£«  T»  3T/3x,  o,^)  ,  (18) 


where  it  is  obvious  that  equations  (14)  through  (18)  satisfy  the  principle 
of  equipresence  [23].  Equations  (10)  through  (12)  and  (14)  through  (18) 
describe  eight  +n  equations  in  the  eight  +n  field  variables  u,  £,  z,  e, 
r,  s,  T,  q,  and  described  in  (1)  through  (9).  These  are  adjoined  with 
boundary  conditions  on  the  surfaces  x  *  0  and  x  *  L  to  prescribe  the  one¬ 
dimensional  field  problem. 

As  detailed  elsewhere  [17,18,20-22],  the  second  law  of  thermodynamics 
[inequality  (13)]  will  constrain  constitutive  assumptions  (14)  through  (18). 
This  is  accomplished  by  defining  the  Helmholtz  free  energy: 

h  2  h(€,  T,  3T/3x,  ^)  -  a  -  Ts  =>  e  =*  h  -f*  Ts  ,  (19) 


which  together  with  the  Clausuis-Duhem  inequality  will  lead  to  the  conclu¬ 


o 


sions  that 


h  =  h(s,  T,  c^) 


(20) 


s  -  3T  **  s(e,  T,  ak) 


3  Z  ,  and 


,  3T  .  n  ,  3T 

q  =  -k  - -  +  0  - 

ox^  cx^ 


where  k  is  the  coefficient  of  heat  conduction  in  the  x^  coordinate  direc¬ 
tion.  Therefore,  equations  (19)  through  (23)  replace  equations  (14)  through 
(18)  as  more  concise  statements  of  the  constitutive  behavior,  and  it  can 
be  seen  that  specification  of  the  Helmholtz  free  energy  will  complete  the 
description  of  the  field  problem. 

Combination  of  equations  (12)  and  (19)  through  (23)  will  result  in  the 


coupled  heat  conduction  equation: 


3h  •  32h  • 

30^  ak  "  '  \ 


-  pT  T-rrsr  £  -  OT 


A  T  +  ^  =  cr  , 
3T  j  x 


where  summation  on  the  range  of  k  is  implied. 

Henceforth  in  this  investigation  it  will  be  assumed  that  there  is 
no  internal  heat  source  (other  than  material  dissipation)  so  that  r  =  0 
in  equation  (24).  In  addition,  it  will  be  assumed  that 
boundary  conditions  are  applied  in  such  a  way 

that  heat  flux  is  negligibly  small  and  q  =  0  in  equation  (24) .  This  last 
assumption  is  not  valid  under  most  physical  circumstances.  However,  it 
can  be  said  that  on  the  basis  of  heat  conduction  equation  (24)  neglecting 
heat  flux  will  result  in  an  upper  bound  for  the  temperature  rise  during 
mechanically  induced  energy  dissipation.  Inclusion  of  this  term  results 
in  a  spacially  dependent  boundary  value  problem  which  is  beyond  the  scope 
of  the  current  research.  However,  the  one-dimensional  model  proposed  herein 


does  encompass  the  heat  flux  phenomenon,  and,  as  such,  will  be  the  subject 
of  a  future  paper  by  the  author. 

Development  of  Governing  Equations  from 
Additional  Constitutive  Assumptions 


In  order  to  construct  the  Helmholtz  free  energy  function  the  elastic 
strain  is  first  defined  to  be 


where  is  the  total  inelastic  strain  in  the  x1  coordinate  system  [24], 

y.  is  the  coefficient  of  thermal  expansion  in  the  x^  coordinate  direction, 

and  9  =  T  -  T  ,  where  T  is  the  initial  temperature  at  which  no  strain  is  ob- 
R  R 

served  under  zero  mechanical  load.  The  inelastic  strain  will  be  discussed 
in  greater  detail  in  the  next  section. 

It  is  now  postulated  that  the  Helmholtz  free  energy  may  be  expanded 
about  the  initial  configuration  in  terms  of  elastic  strain  and  temperature 
as  follows: 


(26) 


where  the  subscript  R  denotes  the  equilibrium  value,  and 

=  free  energy  in  state  R  =  constant,  (27) 

E  =  Young’s  modulus  in  the  x^  coordinate  system,  (28) 

2  2 

=  -T(3  h/3T  )  «  specific  heat  at  constant  volume.  (29) 

E 

Note  that  although  the  first  order  terms  in  t  and  9  have  been  neglected 
due  to  symmetry  conditions  due  to  the  form  of  equation  (25)  coupling  is 
retained  between  total  strain,  inelastic  strain  and  temperature.  Note  also 
that  the  energy  dissipation  due  to  microstructural  change  has  been  neglected 
in  free  energy  equation  (26)  because  this  mechanism  has  been  shown  to  contribute 


/ 


only  a  small  portion  of  energy  (< 10%)  to  the  dissipation  process  [25]. 

Further,  the  fracture  energy  loss  due  to  microvoid  growth,  grain  boundary 
sliding,  and  intergranular  macrofracture  is  neglected  due  to  the  small 
strains  considered  herein. 

Although  the  second  order  Taylor  series  expansion  of  the  Helmholtz 
free  energy  given  in  equation  (26)  may  not  be  adequate  for  characterizing 
the  response  of  many  materials,  it  will  be  shown  in  the  next  section  that 
the  above  equations  are  a  suitable  framework  for  describing  the  material 
behavior  of  the  class  of  materials  considered  herein. 

Substitution  of  equation  (26)  into  energy  balance  law  (24)  and  uti¬ 
lizing  equation  (25)  will  result  in  the  coupled  heat  equation: 

[(Ee  -  Ea  +  EaT_)a,  +  Ea  IT]  -  EaTe  -  cC  T  =  0  ,  (30) 

1  R  1  v 

where  the  terms  in  brackets  arise  due  to  inelastic  response  and  the  following 
term  is  the  classical  elastic  coupling  term  [3].  Equation  (30)  may  be  written 
in  the  following  equivalent  form: 

—  •  —  • 

•  (Ee  -  Ea,  +  EaT  )a,  -  EaTe 


v 


In  order  to  obtain  the  stress-strain  relation  the  Helmholtz  free  energy 
equation  (26)  may  be  substituted  into  equation  (22)  to  obtain 

a  =  E(e  -  a1  -  a9)  .  (32) 

Equations  (31)  and  (32)  together  with  internal  state  variable  growth 
laws  (18)  will  be  sufficient  to  characterize  the  response  of  the  uniaxial 
bar  subjected  to  uniaxial  homogeneous  mechanical  loading  considered  herein. 

Selected  Problems  and  Numerical  Results 
It  has  been  shown  that  stress-strain  relation  (32)  together  with  in- 


ernal  state  variable  growth  laws  (18)  are  equivalent  to  several  models 


recently  proposed  for  thermoviscoplastic  metals  [24].  These  include  Cer¬ 
nocky  and  Krempl  [11,26],  Valanis  [27],  Krieg,  et  al.  [28],  and  Allen  and 
Haisler  [29].  It  can  also  be  shown  that  several  others  are  in  accordance 
with  the  model  developed  herein  [30-34].  To  illustrate  this  point  two 
models  have  been  selected  for  further  discussion. 

Cernockv  and  Krempl 's  stress-strain  relation  may  be  written  in  the 
following  uniaxial  form: 

:  +  K(C,  e,  T)a  -  G(s,  t)  +  m(j,  e,  t) [s  -  ai]  ,  (33) 

where 

E  =  M/K  ,  (34) 

and  parentheses  imply  dependence  on  the  current  values  of  the  quantities 
enclosed.  Equations  (33)  and  (34)  can  be  shown  to  be  in  agreement  with 
stress-strain  equation  (32)  by  defining  the  inelastic  strain  ^  such  that 

-  [?  -  G(e,  T )J/M(a,  s, T )  ,  (35) 

so  that 

t 

(t)  =  /  ot^ (t ' )dt '  ,  (36) 

where  tR  is  the  reference  time,  t'  is  a  dummy  variable  of  integration,  and 
t  is  the  time  of  interest.  Thus,  since  G,  K  and  M  are  not  history  depen¬ 
dent,  Cernocky  and  Krempl' s  model  is  a  single  internal  state  variable  model 
and  equations  (31),  (32),  (35)  and  (36)  describe  the  uniaxial  bar  problem 
using  Cernocky  and  Krempl' s  model. 

To  illustrate  this  point  an  example  problem  is  now  considered.  Sev¬ 
eral  uniaxial  bars  composed  of  mild  steel  are  subjected  to  constant  strain 
rates  at  room  temperature  with  material  properties  as  described  in  Table 
2  of  reference  [11]-  Stress-strain  behavior  and  resulting  temperature 
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rise  are  shown  in  Figs.  2  and  3.  These  results  were  obtained  by  integrating 
equations  (31)  and  (35)  with  a  stable  and  accurate  Euler  forward  integra¬ 
tion  scheme.  Due  to  the  rate  insensitive  nature  of  mild  steel  at  room  tem¬ 
perature,  the  predicted  results  are  identical  for  strain  rates  ranging  from 
0.001  sec  *  to  1.0  sec”*-.  The  negligence  of  heat  flux  over  such  a  wide 
range  of  strain  rates  is  valid  only  under  adiabatic  conditions. 

It  is  significant  to  note  that  the  results  obtained  in  Figs.  2  and 
3  are  identical  to  those  obtained  by  Cernocky  and  Krempl  [11].  This  is 
due  to  the  fact  that  the  assumed  internal  energy  rate  described  by  equa¬ 
tion  (14)  in  reference  [11]  can  be  obtained  in  the  uniaxial  form  by  util¬ 
izing  equations  (19),  (21)  and  (26)  in  this  paper.  Further,  energy  balance 
equation  (55)  in  reference  [11]  can  be  shown  to  be  identical  to  equation 
(26)  derived  herein  by  substituting  equation  (32)  in  this  paper.  Finally, 
it  should  be  pointed  out  that  under  non-adiabatic  conditions  neglecting 
the  heat  flux  in  the  results  obtained  herein  causes  increasing  overesti¬ 
mation  of  the  temperatures  shown  in  Fig.  3  as  the  input  strain  rate  decreases 

Bodner  and  Partom's  model  [35,36]  may  also  be  written  in  the  uniaxial 
form  described  by  equation  (32)  ,  where 


where  Dq  and  n  are  experimentally  obtained  material  constants  and 


/32"Zl\r 

ot,  =  m(Z1  -  a 2)a  -  AZ^I - - - J 


(37) 


(38) 


where  a ^  is  an  internal  state  variable  representing  drag  stress  and  m, 

ZI  ’  Zl*  an<^  r  are  exPetimentally  determined  material  constants.  Although 


equation  (38)  contains  stress  j ,  it  can  be  ’written  in  the  form  described 


ABSOLUTE  STRAIN  (ca/ca> 


(■o/M)  N1VH1S 


by  equation  (18)  by  direct  substitution  of  equation  (32).  Thus,  3odner 
and  Partom's  model  contains  two  internal  state  variables  in  the  form  de¬ 
scribed  above. 

Results  are  shown  in  Figs.  4  through  7  for  uniaxial  bars  of  IN 
100  pulled  at  various  constant  strain  rates  at  an  initial  temperature  of 
1005 °K  (1350°F).  Experimental  data  were  obtained  from  Reference  [37], 
and  the  material  constants  described  above  are  tabulated  in  Reference  [38]. 
The  stress-strain  curves  shown  in  Fig.  4  are  identical  to  those  previously 
obtained  [38]. 

As  a  second  example  using  Bodner  and  Partom's  model  a  uniaxial  bar 
of  IN  100  with  material  parameters  as  described  in  Refs.  [37]  and  [38] 
is  subjected  to  the  cyclic  strain  history  shown  in  Fig.  8  and  at  initial 
temperature  1005°K  (1350°F) .  Analytic  stress-strain  behavior  is  compared 
to  experiment  in  Fig.  8  and  predicted  temperature  change  is  shown  in  Fig.  9 . 

Finally,  a  uniaxial  bar  is  subjected  to  the  multicycle  test  described 
in  Fig.  10,  with  resulting  temperature  rise  shown  in  Fig.  11  It  is  observed 
that  the  model  predicts  a  mean  temperature  rise  of  approximately  3.7°K 
(6.7°F)  per  cycle.  The  linear  increase  in  mean  temperature  with  time  is 
predicted  due  to  the  cyclic  saturation  of  the  material  on  the  second  cycle, 
which  is  in  agreement  with  experimental  observations  at  elevated  temperature. 

Co nclusion 

A  model  has  been  presented  herein  for  predicting  the  temperature  rise 
in  uniaxial  bars  composed  of  thermoviscoplastic  metallic  media.  The  model 
is  also  applicable  to  raultiaxial  conditions,  and  this  has  been  reported  to 
some  extent  in  reference  [20].  Although  the  procedure  used  here  differs 
from  that  proposed  in  reference  [11],  it  has  been  shown  that  the  predicted 
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temperature  change  is  identical  to  results  obtained  by  Cernocky  and  Krempl 
when  their  mechanical  constitutive  equations  are  used.  However,  it  has  been 
shown  herein  that  the  introduction  of  internal  state  variables  leads  to  a 
more  general  model  which  may  be  used  with  virtually  any  thermoviscoplastic 
model  currently  used  for  metals  [24]. 

It  has  been  found  in  the  current  research  that  significant  heating  may 
occur  under  adiabatic  conditions,  expeciallv  during  cyclic  loading,  in  ther¬ 
moviscoplastic  metallic  media.  The  significance  of  this  heating  is  compounded 
by  the  fact  that  material  properties  often  become  extremely  sensitive  in  the 
inelastic  range  of  behavior.  This  issue  has  not  been  considered  herein,  but 
it  certainly  warrants  study  when  transient  temperature  models  become  available. 

Two  important  questions  have  not  been  answered  in  this  research:  1) 
what  effect  does  the  inclusion  of  the  heat  flux  term  have  on  the  predicted 
results;  and  2)  what,  if  anything,  does  the  present  model  have  to  do  with 
experimentally  observed  results?  The  first  question  can  only  be  addressed 
if  spacial  variation  is  admitted  in  the  field  parameters.  The  author  is  cur¬ 
rently  studying  this  question  and  hopes  to  present  results  in  a  future  paper. 
The  second  question  cannot  be  answered  at  this  time  since  it  requires  ex¬ 
tremely  sophisticated  experimentation.  Although  experimental  results  have 
been  obtained  detailing  heat  generation  in  inelastic  media,  it  is  not  pos¬ 
sible  to  compare  the  current  model  since  additional  complex  tests  must  be 
performed  in  order  to  characterize  the  thermoviscoplastic  material  parameters. 
The  author  also  hopes  to  address  this  issue  in  a  future  paper. 
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In  this  paper  a  new  and  efficient  alternative  to  subincrementa- 
tion  is  developed  for  analysis  of  solid  media  with  rate  independent 
elastic-plastic  material  behavior.  This  alternative  method  is  not 
unlike  the  subincrementation  procedure  in  that  it  represents  an  Euler 
integration  of  the  nonlinear  constitutive  equations.  However,  it 
takes  advantage  of  the  fact  that  the  Euler  integration  procedure  as¬ 
sumes  proportional  loading  steps  so  that  when  the  uniaxial  stress- 
strain  curve  is  idealized  as  a  piecewise  linear  relation  very  large 
forward  integration  steps  give  accurate  results.  The  new  procedure, 
which  we  call  the  zeta  method,  is  equally  appropriate  for  cyclic 
loading  with  combined  isotropic  and  kinematic  hardening.  However, 
due  to  the  nonuniqueness  of  the  monotonic  uniaxial  stress-strain  re¬ 
lation  in  rate  dependent  media,  the  method  is  not  appropriate  for  use 
in  viscoplastic  media. 

Although  the  algorithm  deals  only  with  the  evaluation  of  a  clas¬ 
sical  plasticity  based  constitutive  law,  numerical  results  are  re¬ 
ported  herein  for  an  assortment  of  problems  by  the  finite  element 
method.  It  is  shown  via  these  results  that  the  zeta  method  discussed 
herein  provides  not  only  accuracy  which  is  superior  to  the  subincre¬ 
mentation  method,  but  the  resulting  algorithm  also  shows  improved 
numerical  efficiency. 


•n 
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INTRODUCTION 


In  recent  years  the  analysis  of  elastic-plastic  solids  which  behave 
according  to  classical  rate-independent  incremental  plasticity  consti¬ 
tutive  models  has  become  quite  commonplace.  By  far  the  most  often  used 
method  for  nonhomogeneous  boundary  value  problem  solutions  is  the  finite 
element  method.  By  the  nature  of  the  kinematics  and  material  behavior 
this  is  a  nonlinear  field  problem,  and  a  considerable  body  of  research 
has  been  generated  dealing  with  efficient  numerical  solution  of  the 
nonlinear  boundary  value  problem.  One  efficiency  measure  adopted  solely 
to  improve  material  models  and  independent  of  large  deformation  behavior 
is  the  subincrementation  method  [1-3]. 

In  this  paper  the  subincrementation  method  will  be  briefly  reviewed 
and  an  alternative  procedure  will  be  proposed.  It  will  be  shown  that 
this  new  method  gives  both  improved  accuracy  and  efficiency  over  the  sub¬ 
incrementation  method,  and  this  contention  will  be  supported  via  several 
numerical  results. 

REVIEW  OF  THE  ELASTIC-PLASTIC  FIELD  PROBLEM 

An  elastic-plastic  medium  subjected  to  an  isothermal  loading  must 
in  general  satisfy  the  following  conditions  at  all  points  in  its  interior 
V  and  on  its  surface  B:  (1)  conservation  of  linear  and  angular  momentum, 
(2)  conservation  of  mass;  (3)  strain-displacement  relations;  (4)  the 
first  and  second  laws  of  thermodynamics;  and  (5)  stress-strain  relations 
(constitution).  Of  course,  there  are  additional  constitutive  relations 
[4],  but  these  need  not  be  considered  in  order  to  characterize  the  mech¬ 
anical  response  when  internal  heat  generation  is  negligible.  It  will 


I 

be  assumed  chat  Che  above  restriction  holds  for  the  media  considered 
herein,  although  the  method  to  be  considered  here  applies  equally  to  tran¬ 
sient  temperature  phenomena. 

Condition  (5)  is  the  main  topic  of  the  discussion  herein.  In  this 
paper  we  consider  only  elastic -plastic  media:  that  class  of  materials 
for  which  the  stress  (or  strain)  tensor  can  be  considered  to  be  a  time 
independent  functional  of  the  strain  (or  stress)  tensor,  that  is,  stress 
is  dependent  on  the  entire  history  of  strain  but  independent  of  the  time 
scale.  It  has  been  shown  that  this  functional  form  can  often  be  written 
in  an  equivalent  differential  equation  form  [5].  One  common  strain  for¬ 
mulation  is  of  the  general  type  [6] 


da. .  =  C  ds 
ij  ljkl  kl 


(1) 


where  for  small  deformation  0..  and  e. .  represent  infinitesimal  stress 

ij  ij 

and  strain  tensors,  respectively,  and  called  the  effective  mod¬ 

ulus  tensor,  is  the  repository  for  history  dependence  via  its  dependence 


_P 

on  the  equivalent  uniaxial  plastic  strain,  £  ,  which  is  a  metric  in  the 
space  of  plastic  strain  defined  by 

£P(t)  e?.(t) _ 

o  o 

where  t  represents  time  and 


d£.  .  =  d£.  .  -  D.  .dJ, 
lj  ij  ijkl  kl 


(3) 


and  D  ^  is  the  elastic  modulus  tensor,  assumed  to  be  independent  of 

deformation.  In  a  uniaxial  sense  the  effective  modulus  tensor  C..,  , 

ij  kl 


may  be  thought  of  as  a  secant  modulus. 


It  should  be  noted  that  for  finite  deformation  equations  (1)  may 


still  be  applicable  if  .  and  are  replaced  by  frame  indifferent 

ij  iJ 

quantities  consistent  with  hypoelasticity  [7].  Note  also  that  under 
certain  simplifying  assumptions  equations  (1)  reduce  to  the  classical 
Prandtl-  Reuss equations  [8,9]. 

It  has  been  shown  that  equations  (1)  are  consistent  with  the  first 
and  second  laws  of  thermodynamics  [10-12].  Therefore,  since  the  energy 
balance  is  trivially  satisfied  in  an  isothermal  domain  with  negligible 
internal  heat  generation  and  conservation  of  mass  is  satisfied  if  the 
density  is  timewise  constant  during  infinitesimal  deformation,  we  need 
consider  only  conditions  (1) ,  (3)  and  (5)  here. 

In  the  finite  element  method  conditions  (1)  and  (3)  are  usually 
satisfied  via  an  incremental  variational  principle  integrated  over  the 
domain  of  interest.  The  discretization  process  then  entails  reducing 
the  volume  integral  for  the  variational  principle  to  some  sub-domain 
aptly  called  a  finite  element.  The  integration  over  the  volume  of  each 
element  is  usually  sufficiently  difficult  to  require  numerical  integra¬ 
tion,  and  for  this  purpose  a  quadrature  procedure  is  generally  employed 
Thus,  the  integration  is  reduced  to  evaluation  of  the  integrand  at  a 
finite  number  of  integration  points  within  an  element. 

The  volume  integration  of  each  element  yields  a  set  of  matrix  equa¬ 
tions  which  are  assembled  into  a  global  set  of  matrix  equations.  These 
equations  are  nonlinear  even  during  infinitesimal  deformations  due  to 
the  nonlinearity  of  equations  (L).  Therefore,  an  iterative  technique 
is  generally  used  to  solve  the  global  equations  on  each  load  step,  and 
the  constitutive  equations  (1)  must  be  solved  several  times  for  all 


integration  points  on  each  load  step  until  convergence  occurs. * 


SOLUTION  OF  THE  CONSTITUTIVE  EQUATIONS 


It  has  become  common  practice  in  the  literature  to  incrementalize 
equations  (1)  by  simply  replacing  differentials  d~_  and  d£„  with  in¬ 
crements  Aj^.  and  ,  respectively.  However,  since  the  integrand 

depends  on  the  strain  tensor  during  the  load  increment,  an  important 
task  becomes  the  integration  of  equations  (1)  over  some  input  increment 
in  the  strain  tensor,  viz.: 

a  . . (t+At) )  £, . (t+At) 

r  r1: 


\  ~  _ 

.  .  — 
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C.  .  (£P)d£, 

ijkl 


kl 


where  t  represents  the  time  at  the  start  of  a  load  increment,  and  t+At 
is  the  time  at  the  end  of  a  load  increment. 

Equations  (4)  definitely  present  a  uniqueness  problem  since  the 
strain  tensor  may  be  cycled  during  the  time  increment  At.  In  order  to 
avoid  this  difficulty  a  sufficient  condition  may  be  adopted  which  is 
not  unlike  the  condition  required  to  obtain  the  Mises-Hencky  deforma¬ 
tion  theory  [15]  from  the  Prandtl-Reuss  equations.  It  is  assumed  that 
during  the  time  increment  At  all  components  of  the  strain  tensor  increase 
monotonically  via  the  relation 

— P 

dc^  =  K^de  »  =  constants  ,  (5) 

— P 

where  d£  must  be  a  monotonically  increasing  function  of  strain  during 
plastic  loading  over  the  time  increment  At.  Substitution  of  equations 
(5)  into  (4)  yields 

*For  a  more  complete  discussion  of  the  finite  element  method  applied  to 
elastic-plastic  media,  see  references  [13]  and  [14]. 
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which  is  obviously  a  unique  relation  due  to  the  monotinicity  of  the 
plastic  strain  increment  during  the  time  increment  -It.  Now  define  the 
following  fourth  order  tensor  which  is  constant  over  the  time  step  it: 


^(t+it) 


r ' 

ljkl 


-I*  /  ^ 

_ p  ** 

e  (t) 


(rp)dip 


(7) 


Then  substitution  of  definitions  (7)  into  equation  (6)  gives 


10..  =  C !  . .  ,  K.  . 
ij  ljkl  tel 


=  C! 


(8) 


ijkl  kl 

which  is  the  exact  incremental  relation  which  should  be  used  with  the 

incremental  variational  principle.  It  should  be  noted  that  equations 

(1)  and  (8)  are  by  no  means  equivalent  since  C.'.,  ,  can  be  seen  from 

ij  kl 

definitions  (7)  to  represent  an  average  effective  modulus  tensor  during 
the  time  increment  it.  Unfortunately,  equations  (7)  cannot  be  integrated 

— p 

precisely  because  the  upper  limit  of  integration  £  (t+it)  cannot  be  de¬ 
termined  until  equations  (8)  have  been  evaluated. 

It  will  be  recalled  that  the  equivalent  uniaxial  plastic  strain 
can  be  shown  to  be  equal  to  the  axial  plastic  strain  when  a  bar  is  pulled 
uniaxially  [16].  Now  define  the  equivalent  uniaxial  stress 

J. 


f3  ,  KK  .  w  ^KK  -  ,T* 

(.7  r:'J  ■ 


(9) 


which  can  also  be  seen  to  be  equivalent  to  the  axial  stress  when  a  bar 
is  pulled  uniaxially  [16],  Thus,  the  information  required  to  charac¬ 
terize  equations  (6)  is  obtainable  from  a  single  monotonically  increasing 
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equivalent  uniaxial  plastic  strain  diagram  as  shown  in  Fig.  1.  N'ote 
that  the  curve  shown  represents  nothing  more  than  a  uniaxial  stress- 
strain  diagram  with  abscissa  transformed  via  definition  (3).  In  addi¬ 
tion,  for  combined  isotropic-kinematic  hardening  the  ordinate  in  Fig. 

1  should  be  transformed  as  well  [17]. 

It  is  apparent  from  Fig.  1  that  for  a  continuously  work  hardening 
material  the  relation  between  the  equivalent  uniaxial  plastic  strain 

_ p 

£  and  the  slope  of  the  uniaxial  stress-plastic  diagram  is  unique,  that 


is, 


-rP 


(10) 


where  F  is  a  bijective  mapping.  Therefore,  the  effective  modulus  tensor 
may  be  written  alternatively  as 


(ID 


Thus,  because  the  effective  modulus  tensor  C...  ,  is  a  nonlinear  function 

ljkl 

_  — P 

of  the  plastic  secant  modulus  d<J/d£  ,  integration  of  equations  (7)  is 
not  a  trivial  task.  In  order  to  avoid  this  numerical  complexity  it  is 
not  uncommon  to  simply  approximate  the  effective  modulus  tensor  by 


Cijkl(e  }  ~  Cijkl(e  (C)) 
thus  reducing  integration  of  equations  (6)  to 

-hj  a  ■  cijki(=P<t))“ 

which  is  Euler's  method  of  forward  integration, 
is  nothing  more  than  a  simple  first  order  Taylor 


kl 

Obviously,  since 
series  expansion 


(12) 


(13) 


this 

its 


accuracy  will  depend  on  the  relative  nonlinearity  or  curvature  of  the 
uniaxial  stress-plastic  strain  curve  during  a  given  load  increment* 
This  condition  is  illustrated  for  the  uniaxial  case  in  Fig.  2. 


In  many  computer  codes  the  uniaxial  stress-strain  diagram  is  input 
in  piecewise  linear  fashion  as  shown  in  Fig.  3.  While  it  is  not  clear 
that  this  piecewise  linearization  is  motivated  by  anything  beyond  sim¬ 
plification  of  input  data,  it  has  the  added  benefit  that  it  helps  alle¬ 
viate  the  numerical  integration  problem  noted  above  and  described  in 

— p 

Fig.  4.  In  fact,  so  long  as  the  plastic  strain  increment  Lz  does  not 
subtend  a  slope  discontinuity  during  a  load  step  equations  (11)  indi¬ 
cate  that  approximation  (12)  will  reproduce  the  piecewise  linear  curve 
exactly.  Therefore,  the  accuracy  of  equations  (13)  is  limited  only  by 
the  accuracy  with  which  one  can  reproduce  the  exact  curve  of  stress  versus 
plastic  strain  with  a  piecewise  linear  curve.  Mathematically,  the  slope 
continuity  condition  is  satisfied  if  one  can  find  values  of  the  equiva- 

— p 

lent  plastic  strain  at  slope  discontinuities  ,  as  shown  in  Fig.  3, 
such  that  for  the  current  load  increment 


(14) 


for  all  equivalent  plastic  strains  in  the  range 

~P(t)  <  £P  IP(t+At)  .  (15) 

However,  condition  (14)  cannot  be  a  priori  guaranteed  in  practice  be¬ 
cause  in  a  non-homo geneous  boundary  value  problem  the  equivalent  uniaxial 
plastic  strain  varies  spatially.  Whereas  one  integration  point  may 
undergo  a  very  small  or  even  zero  (elastic)  plastic  strain  increment 
during  a  specified  increment  in  surface  tractions,  another  point  under 
high  stress  concentration  may  undergo  a  plastic  strain  increment  which 
subtends  one  or  more  discontinuities  in  the  piecewise  linear  equivalent 
uniaxial  stress  versus  equivalent  uniaxial  plastic  strain  diagram.  Thus, 
as  can  be  seen  from  definitions  (2)  and  (3)  the  equivalent  uniaxial 
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Fig.  3.  Piecewise  Linearization  of  the  Equivalent 
Uniaxial  Stress  Versus  Equivalent  Uniaxial  Plastic 
Strain  Diagram 


plastic  strain  increment  can  be  determined  only  after  the  stress  incre¬ 
ment  tensor  has  been  found. 

REVIEW  OF  THE  SUB INCREMENT  AT ION  METHOD 


In  order  to  improve  the  accuracy  of  approximation  (13)  subincremen¬ 
tation  has  been  proposed  [1,2].  In  this  method  Stricklin,  et  al.  define 
the  equivalent  uniaxial  strain  increment 


It  =  J  |  dt..d£.. 

>3  IT  ij 


_  (16) 

i-l  iJ 

This  quantity  is  evaluated  over  a  specified  load  step  and  is  then 
compared  to  an  input  parameter  called  the  allowable  total  strain  incre¬ 


ment  (dc^)  as  follows 


M  = 


(17) 


dc 


AL 


where  M  is  rounded  off  to  the  nearest  integer  greater  than  Eero.  Equa¬ 
tions  (6)  are  then  evaluated  M  times  for  the  strain  subincrement 

„  Ac.  . 

Ac .  .  =  -rr1  ,  (18) 

ij  M 

and  on  each  subincrement  the  effective  modulus  tensor  is  updated  to 
reflect  the  current  equivalent  uniaxial  plastic  strain.  Based  on  numer¬ 
ical  evidence,  Stricklin  suggests  that  dc  should  be  no  greater  than 
0.0005  in. /in.,  although  our  experience  indicates  that  values  as  small 
as  .00005  in. /in.  are  sometimes  required  to  maintain  accuracy  of  solu¬ 
tion. 


In  order  to  illustrate  the  effect  of  subincrementation  let  us  ex¬ 
amine  a  single  example.  Suppose  we  consider  a  bar  subjected  to  a  grad¬ 
ually  increasing  homogeneous  uniaxial  stress  state.  Because  conditions 


13  - 


(1)  and  (3)  are  satisfied  trivially  we  need  only  consider  approximate 
constitutive  equations  (13).  Since  the  input  material  properties  will 
be  described  via  a  piecewise  linear  equivalent  uniaxial  stress  versus 
equivalent  uniaxial  plastic  strain  diagram  as  shown  in  Fig.  3,  and  be¬ 
cause  this  boundary  value  problem  is  equivalent  to  the  experiment  which 
produced  the  material  input  data,  an  exact  analysis  using  equations 
(13)  should  reproduce  Fig.  3  precisely.  In  fact,  using  subincrementa¬ 
tion  will  yield  the  results  shown  in  Fig.  4  when  a  single  slope  discon¬ 
tinuity  is  encountered  in  a  given  load  step.  It  can  be  seen  from  the 
figure  that  the  total  error  is  incurred  during  the  plastic  strain  sub¬ 
increment  subtending  the  slope  discontinuity.  The  effect  then  of  sub¬ 
incrementation  is  simply  to  improve  the  approximate  integration  of  equa 
tions  (6). 

A  PROPOSED  MODIFICATION 

It  is  apparent  from  the  above  discussion  that  subincrementation 
will  often  require  multiple  evaluations  of  equations  (13)  for  each  inte 
gration  point.  Since  these  equations  must  be  evaluated  at  each  inte¬ 
gration  point  in  the  body  and  often  several  times  for  each  load  step 
in  order  to  obtain  equilibrium  convergence,  considerable  computational 
time  can  be  spent  in  this  process.  Detailed  herein  is  a  numerical  pro¬ 
cedure  for  integrating  equations  (7)  which  is  both  more  accurate  and 
more  computationally  efficient  than  subincrementation.  We  call  this 
method  the  zeta  method. 

The  method  proposed  here  is  a  simple  extension  of  a  procedure  uti¬ 
lized  by  Krieg  and  Duffey  [18]  for  the  transition  step  from  elastic  to 
elastic-plastic  material  behavior.  The  primary  extension  is  that  each 


subsequent  slope  discontinuity  in  the  piecewise  linear  equivalent  uni¬ 
stress  versus  equivalent  uniaxial  plastic  strain  diagram  is  treated  ex¬ 
actly  like  a  subsequent  yield  surface.  Since  equations  (13)  are  exact 
under  conditions  (14)  and  (13) ,  no  subincrementation  will  be  required 
to  obtain  precise  results  between  slope  discontinuities. 

In  order  to  see  how  the  zeta  method  works,  consider  a  material  point 
which  is  in  a  post-yielded  state  at  time  t  and  with  equivalent  uniaxial 

— p 

plast  strain  e  (t) ,  as  shown  in  Fig.  o. 

According  to  Krieg  and  Duffey  [18],  the  value  of  the  stress  tensor 
at  the  material  point  necessary  to  bring  the  equivalent  uniaxial  stress 
state  to  the  i+lth  slope  discontinuity  is  defined  by 

j1*1  =  0.  1  +  ;  Ac  (19) 

ij  ij 

where  Ac„  the  increment  predicted  by  equations  (13)  using  the  input 

total  strain  increment  As..  and  (  is  a  scalar  factor  to  be  determined. 

ij 

In  order  to  determine  zeta  definitions  (19)  are  substituted  into 
the  yield  criterion  used  in  the  model.  For  example,  if  von  Mises' 
yield  criterion  is  utilized,  equations  (19)  will  result  in 


4  a . . t  +  ;Ac . ,  - 
2  JJ  J-j 


a. . 
J  ij 

J  W 


t  (aKK+CA<JKK)  - 

t  j  L.J  .  .  -  - r -  C  .  . 

ij  3  ij 


<=W‘!-  <20> 


Solving  the  above  equation  for  zeta  will  result  in 


r  -B  +  V(B  -  4 AC) 
2A 

where 


A  =  Ac . . Ac . .  , 

ij  ij 


ij  ij 
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and 


C  *  j.  .a.  .  -  (a  ,  )**  .  (24) 

ij  il  1+1 

Utilizing  equations  (21)  through  (24)  the  value  of  zeta  may  be 
obtained.  If  zeta  is  greater  than  or  equal  to  unity,  the  input  total 
strain  increment  will  not  take  the  equivalent  uniaxial  stress  beyond 
the  next  slope  discontinuity  and  the  results  of  equations  (13)  may 
be  considered  correct.  In  other  words 

r  >  1  =>  Aa . .  =  Ac. .  .  (25) 

“  ij  iJ 

If,  on  the  other  hand,  zeta  is  less  than  unity,  the  predicted  stress 
increment  tensor  is  incorrect  and  equations  (13)  must  be  modified. 

This  is  accomplished  by  first  constructing  the  input  strain  increment 
necessary  to  bring  the  equivalent  uniaxial  stress  to  the  slope  discon¬ 
tinuity: 

Asi+1  -  5  Ae. .  ,  (26) 

ij  il 

where  As.  .  is  the  input  strain  increment.  The  values  of  As^T  are 
ij  iJ 

then  substituted  into  equations  (13)  to  produce 

"i?  -  •  (27) 

the  remaining  portion  of  the  stress  increment  tensor  is  calculated 
by  first  determining  the  remainder  of  the  input  strain  increment 

As*.  =  (1  -  5)Aeij  ,  (28) 

noting  from  definitions  (26)  and  (28)  that 


As. .  =  Asi+1  +  As*.  .  (29) 

ij  ij  ij 

£ 

The  remainder  of  the  strain  increment  tensor  As^  is  then  substituted 
into  equations  (13)  to  give  the  remainder  of  the  stress  increment  tensor 


17  - 


AcR.  =  c  (cP  )AcR. 
ij  ijkl  i+l  ij 


(30 


Thus,  the  total  stress  increment  tensor  is  given  by 


Ac..  =  Ac 1+1  +  AcR 
13  ij  tJ 


(31 


It  is  easily  verified  that  the  above  procedure  will  result  in  an  equiv¬ 
alent  uniaxial  stress  and  equivalent  uniaxial  plastic  strain  [utilizing 
definitions  (2)  and  (3)]  which  lie  on  the  equivalent  uniaxial  stress 
versus  equivalent  uniaxial  plastic  strain  diagram.  It  should  also  be 
pointed  out  that  although  the  actual  yield  surface  is  updated  through¬ 
out  plastic  loading,  the  equivalent  uniaxial  stresses  corresponding  to 
slope  discontinuities  should  at  no  time  be  altered. 

Although  the  above  procedure  has  been  discussed  here  only  in  the 
context  of  isotropic  hardening,  it  is  also  applicable  to  more  complex 
yield  criteria  and  work  hardening  rules  [19,20], 


COMPUTER  CODE  FLOWCHART 


The  following  chart  outlines  in  abbreviated  form  the  application 

of  the  c  method  for  a  given  increment  in  the  total  strain  tensor  Ae.. 

ij 

— P 

and  equivalent  uniaxial  plastic  strain  at  the  start  of  the  step  £  (t). 

a)  Set  Ac. .  =  0 . 

ij 

— p 

b)  Evaluate  C.  .(e  (t)). 

IJ  K.1 

c)  Obtain  Ac„  using  equations  (13). 

d)  Determine  c.  ,  from  Fig.  5. 

l+l 

e)  Calculate  A,  B,  and  C  using  equations  (22)  through  (24). 

f)  Determine  zeta  using  equation  (21). 


g)  If  (  A  1  go  to  step  q) . 


h)  If  £  <  1  evaluate  A;  _  using  equations  (26). 

i)  Determine  Aj^T^  using  equations  (27). 

j )  Set  Ac . .  =  Ac . .  +  Aji+1 . 

ij  ij  ij 

k)  Calculate  At^  using  equations  (28). 

l)  Set  Ac . .  =  Ac . . . 

ij  ij 

P 

m)  Calculate  Ae„  using  equations  (3). 

— p 

n)  Determine  Ac  using  equation  (2). 

— P  _p  _p 

o)  Set  c  (t)  =  c  (t)  +  Ac. 

p)  Go  to  step  b) . 


1)  Set  Ac. .  =  Ac. .  +  Ac. . . 

ij  iJ  iJ 


DISCUSSION  OF  RESULTS 

In  this  section  the  results  generated  using  the  zeta  method  as 
well  as  the  subincrementation  method  for  solving  the  constitutive  equa¬ 
tions  of  classical  plasticity  are  presented.  Both  solution  techniques 
have  been  incorporated  into  a  finite  element  computer  program  which 
uses  constant  strain  triangular  elements.  The  formulations  have  been 
cast  into  a  2-dimensional  plane  stress  format. 

In  order  to  compare  the  efficiency  of  the  two  different  methods 
their  respective  solution  times  will  be  compared.  This  was  accomp¬ 
lished  by  using  the  built  in  timer  (clock)  subroutine  used  in  the 
Fortran-H  Extended  language  available  on  the  AMDAHL  470  V6  located 
on  the  Texas  A&M  University  campus.  Two  times  will  be  given  in  the 
analysis:  1)  time  spent  using  the  constitutive  package,  and  2)  total 

time  spent  in  solving  the  specified  boundary  value  problem. 

Two  boundary  value  problems  have  been  selected  for  comparing  the 
two  constitutive  packages:  a  highly  yielded  uniaxial  bar  subjected 


•s. 


^ ' 


Wv*1 
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Co  uniaxial  tension  only  and  a  thick  walled  pressure  vessel  subjected 
to  internal  pressure  sufficient  to  yield  a  broad  band  of  the  vessel. 
Both  specimens  are  assumed  to  be  made  of  5Q86-H34  aluminum  with  piece- 
wise  linearized  room  temperature  properties  shown  in  Figure  6.  The 
finite  element  mesh  used  for  the  uniaxial  bar  as  well  as  the  load 
input  diagram  are  shown  in  Figure  7.  Only  two  elements  are  necessary 
to  represent  the  bar  because  the  boundary  value  problem  is  homogeneous. 
However,  if  the  problem  was  inhomogeneous  then  mesh  refinement  would 
be  necessary  to  increase  the  accuracy  of  the  solution.  It  should 
nevertheless  be  pointed  out  here  that  the  number  of  elements  used  in 
the  mesh  is  directly  proportional  to  the  computational  time  required 
in  the  constitutive  package.  Another  factor  influencing  the  computa¬ 
tion  time  is  the  non-linearity  of  the  given  stress  strain  curve.  In¬ 
creasing  the  number  of  piecewise  linearities  in  the  idealized  uniaxial 
stress  strain  curve  will  increase  the  computational  time  required  by 
the  zeta  method.  Although  this  increase  in  piecewise  linearizations 
will  not  greatly  affect  the  computational  time  required  by  the  subin¬ 
crementation,  it  will  have  an  adverse  affect  on  the  accuracy  of  this 
method. 

The  results  of  the  uniaxial  bar  test  are  shown  in  Fig.  8  and 'the 
comparative  solution  times  are  given  in  Tables  1  and  2.  Fig.  3  shows 
the  output  axial  displacement  versus  time  for  the  zeta  method  as  well 
as  the  subincrementation  with  various  allowable  errors  in  equivalent 
uniaxial  strain  shown  in  parentheses.  Table  1  shows  a  comparison  of 
solution  times  for  each  load  step,  while  Table  2  gives  a  more  detailed 
comparison  of  solution  times  for  the  final  time  step.  Several  dif¬ 
ferent  cases  were  run  using  the  subincrementation  code  in  order  to 
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UNIAXIAL  STRESS  STRAIN  CURVE 

5086-H34  ALUMINUM 


STRAIN  IN/IN 


2 


3 


UNIT  TIME 


Lai  Bar  Finite  Element  Mesh 


I 


TABLE  2: 


COMPARISON  OF  SOLUTION  TIMES  FOR 
THE  LAST  LOAD  STEP  ON  THE 
UNIAXIAL  BAR 


SOLUTION  TIMES 

(sec) 

TIME 

ZETA 

METHOD 

SUBINC 

(0.005) 

SUBINC 

(0.0005) 

SUBINC 

(0.0001) 

SUBINC 

(0.00005) 

CONSTITUTIVE 

0.063492 

0.0027 

0.0555 

0. 18764 

0. 37684 

SUBINC/. 063492 

.043 

.716 

2.955 

5.  935 

TOTAL  TIME  (SEC) 

0.0968 

0.00798 

0.07761 

0.22315 

0.40994 

SUBINC/.  0968 

.0824 

.802 

2.305 

4.235 

illustrate  the  difference  in  solutions  and  comparative  times  for  dif¬ 
ferent  values  of  the  allowable  uniaxial  stress  increment  d£  used 
in  determining  the  number  of  subincrements.  In  general,  increasing 
the  required  accuracy  of  the  subincrementation  method  also  increases 
the  computational  time. 

The  results  of  this  test  show  that  in  order  to  obtain  solution 

accuracy  by  subincrementation  which  is  comparable  to  the  zeta  method, 

the  allowable  strain  increment  must  be  of  the  order  of  0.00005  in/in. 

In  fact,  the  larger  the  allowable  strain  increment  the  less  solution 

time  required,  and  in  fact  for  d£  =  0.005  and  0.0005  the  solution 

AL 

time  was  less  for  the  subincrementation  method.  However,  the  resulting 
accuracy  was  very  poor.  Table  1  indicates  that  subincrementation 
(£  =  0.00005  in. /in.)  requires  3.171  times  as  much  computer  time  as 

the  zeta  method  for  the  uniaxial  bar  problem.  Although  this  is  a 
rather  large  difference  in  relative  times,  since  only  a  two  element 
problem  has  been  run,  the  difference  in  actual  cost  is  small.  However, 
on  an  extremely  large  scale  problem  obvious  savings  would  result. 

The  finite  element  mesh  used  for  the  thickwall  pressure  vessel 
is  shown  in  Fig.  9  and  the  load  input  diagram  is  shown  in  Fig.  10. 

The  results  of  this  test  are  shown  in  Fig.  11  as  well  as  Tables  3 'and 
4.  These  results  are  for  the  final  pressure  of  P  =  30,000  PSI.  It 
should  be  noted  that  if  one  applies  increasing  pressure  to  the  speci¬ 
men  then  more  of  the  elements  will  yield  in  the  outer  regions  of  the 
thickwalled  pressure  vessel,  resulting  in  a  higher  solution  time  be¬ 
cause  more  time  is  spent  in  the  constitutive  package  and  more  time 
in  iterating  on  the  correct  nonlinear  solution. 


Fig.  9.  Finite  Element  Mesh  For 
Thick  Walled  Pressure  Vessel 


THICK  WALLED  PRESSURE  VESSEI 


LOCATION  IN 


TABLE  3:  RADIAL  DISPLACEMENTS  AND  SOLUTION 

•  *  « 

TIMES  FOR  THICK  WALLED  PRESSURE  VESSEL 

RADIAL 

LOCATION 

RADIAL 

DISPLACEMENT 

(IN) 

Li) 

s 

'  •v 

ZETA 

METHOD 

SUBINC 

(0.001) 

SUBINC 

(0.0005) 

SUBINC 

(0.0001) 

SUBINC 

(0.00005) 

^■1 

2.0 

0.01324344 

0.01285087 

0.01289475 

0.01292329 

0.0129211 

M 

2.2222 

0.01167517 

0.01133769 

0.01137959 

0.01141027 

0.01140665 

2. 4444 

0.0105329 

0.01026293 

0.01029393 

0.01032073 

0  01031638 

U 

2.6666 

0.00971329 

0.00947785 

0.00960516 

0.00953320 

0.00952741 

3.0 

0.00883162 

0.00862648 

0.00866162 

0.00868949 

0.00868414 

i1' 

•Vrf 

3.  3333 

0.00821108 

0.00802864 

0.00806262 

0.00809172 

0.00808688 

4.0 

0.00728626 

0.00714737 

0.00717933 

0.00720508 

0.00720022 

■■  > 

.<■ 

CONSTITUTIVE 
PACKAGE  TIME  (SEC) 

0.691703 

0.6132077 

0.6854118 

1.056744 

1.503111 

'v 

^  ‘ 

*4. 

y 

SUBINC/. 691703 

.886 

.991 

1.527 

2.  173 

TOTAL  TIME  (SEC) 

3.  1546827 

3.0397627 

3.  1661748 

3.61964 

4.0580 

£ 

SUBINC/3. 1546827 

.  964 

1.003 

1.147 

1.286 

gw 

TABLE  4: 


COMPARISON  OF  SOLUTION  TIMES  FOR 
THE  LAST  LOAD  STEP  ON  THE 
THICK  WALLED  PRESSURE  VESSEL 


SOLUTION  TIMES  (sec) 


TIMES 

ZETA 

METHOD 

SUBINC 

(0.001) 

SUBINC 

(0.0005) 

SUBINC 

(0.0001) 

SUBINC 

(0.00005) 

CONSTITUTIVE 

PACKAGE  TIME  (SEC) 

0.476943 

0. 398788 

0.468156 

0.838656 

1.278211 

SUBINC/.  476943 

.830 

.  982 

1.758 

2.680 

TOTAL  TIME  (SEC) 

1.662752 

1.565564 

1.665976 

2.  143309 

2.579693 

SUBINC/1. 662752 

.941 

1.002 

1.289 

1 . 551 
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The  results  of  the  test  are  basically  the  same  as  those  for  the 
uniaxial  bar.  From  Table  3  it  can  be  seen  that  for  comparable  solu¬ 
tion  accuracy  the  subincrementation  method  takes  1.527  to  2.173  times 

— P 

more  constitutive  time  (depending  on  £  )  than  the  zeta  method  and 

ivL* 

1.147  to  1.286  times  greater  total  time.  Closer  examination  shows 

— p 

that  even  for  the  smallest  £..  =  0.00005  the  solutions  still  differ 

AL 

from  the  zeta  method  and  that  there  is  no  noticeable  difference  be¬ 
tween  solutions  for  dc^T  =  0.0001  and  de^T  =  0.00005.  In  fact,  the 

AL  AL 

— p 

results  tend  to  be  less  accurate  for  dj:  =  0.0000O.  This  can  be  at- 

AL 

tributed  to  numerical  roundoff  error  because  the  increments  in  the 
strain  are  so  small  that  further  improved  accuracy  is  not  possible. 

By  constrast,  there  is  no  numerical  roundoff  error  apparent  in  the 
zeta  method. 

CONCLUSION 


The  objective  of  this  research  has  been  to  produce  an  alterna¬ 
tive  to  subincrementation  which  results  in  a  substantial  improvement 
in  computational  efficiency.  This  new  method  has  been  shown  by  ex¬ 
ample  to  give  not  only  improved  efficiency,  but  also  slightly  greater 
accuracy  of  results.  The  following  general  conclusions  can  be  made: 

1)  in  order  to  produce  results  by  the  subincrementation  method 
which  are  comparable  in  accuracy  to  the  zeta  method,  signi¬ 
ficantly  greater  computation  time  is  required  by  the  former 
me  t  hod , 

2)  increasing  required  accuracy  in  the  allowable  equivalent  uni¬ 
axial  strain  increment  can  lead  to  roundoff  error  when 


subincrementation  is  utilized; 


3)  sub increment a cion  necessarily  produces  errors  in  predicted 
stresses  whenever  a  slope  discontinuitv  is  subtended  in  the 
uniaxial  stress-strain  curve, 

4)  the  zeta  method  follows  the  uniaxial  stress-strain  curve  ex¬ 
act  ly; 

5)  both  subincrementation  and  the  zeta  method  approximate  the 
load  path  to  be  radial  during  each  load  increment, 

6)  piecewise  linearization,  although  merelv  a  numerical  conven¬ 
ience  in  subincrementation,  is  necessary  in  order  to  utilize 
the  zeta  method, 

7)  the  zeta  method  can  be  used  with  cyclic  hardening  models  of 
plasticity,  and 

8)  the  zeta  method  may  not  be  appropriate  for  use  in  rate  depen¬ 
dent  viscoplastic  media. 
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ABSTRACT 


In  recent  years  research  on  the  subject  of  large  space  structures  (LSS)  has 
dramatically  increased.  This  research  is  of  an  extremely  broad  nature,  extending 
from  materials  characterization  to  dynamic  response  as  well  as  active  control. 
This  paper  is  intended  to  be  a  review  of  recent  and  important  advances  in  this 
field  of  research.  The  subject  matter  is  divided  into  six  major  headings:  1) 
expected  applications  of  LSS;  2)  the  environment;  3)  materials;  A)  solution  tech¬ 
niques  for  LSS;  5)  damping  technology;  and  6)  preliminary  design/analysis/experi¬ 
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I.  Expected  Applications  of  LSS 

With  the  STS  operational,  the  opportunity  exists  to  exploit  space  through 
a  myriad  of  enterprises  that  will  be  economically  beneficial  to  mankind.  Studies 
conducted  to  date  point  out  many  projects  which  will  require  the  construction 
of  very  large  structures  in  space.  Possible  applications  of  LSS  include  solar 
energy  collection  and  transmission,  communications,  navigation,  resources  and 
environmental  survey,  border  and  coastal-water  survery,  area  illumination, 
manufacturing,  radio  astronomy,  and  various  scientific  opportunities  that  a 
space  station  will  provide. 

Hedgepeth  [1]  conducted  a  survey  of  future  requirements  for  large  space 
structures  in  1976.  In  it  he  states  that  the  applications  for  LSS  can  be  grouped 
into  three  categories:  (1)  large  surfaces  for  power,  (2)  large  antenna  to  re¬ 
ceive  and  transmit  energy  over  the  rf  loadwidth,  and  (3)  space  platforms  to 
provide  area  for  general  utlizations.  An  example  of  an  application-type  paper 
Woodcock  [2]  would  fall  into  Hedgepeth's  first  category.  Folder  and  Dienemann 
[3],  and  Wannlund  [4]  project  the  need  for  large  antennas.  Holloway  and  Garrett 
[5],  and  Stone  [6]  examine  the  utilization  of  space  platforms.  Hagler  [7] 
addresses  the  feasibility  of  two  possible  applications:  (1)  a  200m  radio 
astronomy  telescope,  and  (2)  a  5000  MW  sps. 

Cord,  Kruszewski,  and  Gustaferro  [8],  and  Daros,  Freitog,  and  Kline  [9] 
project  probable  chronologies  of  development  of  LSS. 

Hedgepeth,  Mikulas,  and  Macneal  [10]  concern  the  economic  feasibility  of 
building  large  structures  in  space. 

Since  1978,  NASA-LaRC  has  sponsored  an  annual  conference  on  LSST.  Papers 
by  James  [11]  and  Cord  [12]  are  introductions  to  the  two  volumes  of  the  1980 
conference  proceedings. 
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II.  The  Environment 


To  be  feasible,  LSS  will  have  a  function  safely,  reliably,  and  effic¬ 
iently  through  several  years  of  service.  They  must  therefore  be  designed  with 
particular  emphasis  given  the  environment  in  which  they  will  be  expected  to 
function.  Preliminary  studies  indicated  that  most,  if  not  all,  applications 
of  LSS  will  place  them  in  GEO.  They  must  therefore  be  designed  for  the  geosyn¬ 
chronous  radiation  environment  and  severe  thermal  cycling. 

Ambient  radiation  is  of  particular  concern  to  the  designer  because  it 
causes:  (1)  a  serious  hazard  for  humans,  (2)  degrading  and  damaging  effects 

on  materials,  (3)  saturation  of  instruments,  and  (4)  interference  with  scien¬ 
tific  measurements  and  observations.  Stassenopoulous  [13]  provides  a  study 
of  the  geostationary  radiation  environment. 

Because  spacecraft  at  GEO  can  charge  to  levels  significant  enough  to  im¬ 
pact  the  degrade  system  performance,  there  is  a  need  for  space  vehicle  design 
engineers  to  have  a  "worst  case"  charging  environment.  Gussenhoven  and  Mullen 
[14]  discuss  an  environment  that  produced  one  of  the  largest  satellite  frame 
potentials  on  the  SCATHA  satellite  in  its  first  year  of  operation.  Presented 
is  an  extremely  severe  geostationary  environment  having  an  exceptionally 
dense  plasma  population  in  the  high  energy  range  (20  keV) .  This  is  the 
only  paper  included  in  this  survey  regarding  spacecraft  charging,  but  much 
has  been  written  and  this  paper  lists  many  references. 

The  success  of  nearly  every  mission  envisioned  for  LSS  is  critically 
dependent  upon  the  maintenance  of  a  stable,  close-tolerance  geomtry.  Thus 
the  prediction  and  control  of  structural  deformations  under  the  influence 
of  the  orbital  environment  become  major  design  requirements  and  the  candi¬ 
dates  structure's  characteristic  responses  become  important  considerations 
in  concept  selection.  Brogren.  Barclay,  and  Straayer  [15]  provide  the  de- 
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signer  with  a  tool  for  making  rapid  estimates  of  the  response  of  these 
structures  to  the  thermal  environment  encountered  in  Earth  orbit.  Thermal 


analysis  of  relatively  sparse  structures  in  the  space  environment  has  cus¬ 
tomarily  omitted  consideration  of  shadowing  by  up-sun  structural  members. 

This  convention  has  been  frequently  questioned  in  the  case  of  lattice-type 
structures  supporting  very  large,  near  planar.  Earth-facing  surfaces.  For 
these,  significant  shadowing  can  occur  whenever  the  solar  vector  is  nearly 
tangent  to  the  orbital  path.  Oneil,  [16]  deals  with  modelling  this  problem 
area . 

Bosma  and  Levadou  [17]  examine  passive  thermal  control  coatings  for 
Spacelab.  This  may  indirectly  be  of  interest  to  those  studying  thermal  control 
for  LSS . 

III.  Materials 

The  materials  used  for  the  construction  of  LSS  will  be  exposed  to  an 
environment  that  will  be  seriously  conducive  to  material  degradation.  The 
geosynchronous  environment  has  an  exceptionally  dense  plasma  population  in 
the  high  energy  range,  and  there  is  little  natural  protection  from  the  solar 
electromagnetic  spectrum.  The  structure  will  be  exposed  to  severe  thermal 
cycling  (one  study  indicates  a  range  of  from  88K  to  377K,  another  from  77K 
to  394K) . 

Rosen  and  Hashin  [18]  have  studied  the  relationship  between  the  consti¬ 
tuent  and  composite  thermoelastic  properties. 

Schapery  [19]  has  determined  bounds  on  coefficients  of  thermal  expansion 
for  materials  based  on  principles  of  thermoelasticity.  This  work  is 


applicable  for  both  istropic  and  anisotropic  composite  materials. 


Halpin  [20]  has  estimated  expansion  and  stiffness  for  oriented  short 
fiber  composites  with  the  use  of  Halpin-Tsai  equations. 

The  effects  of  thermal  expansion  in  metal-matrix  composites  have  also 
been  investigated  by  Wang  and  Sutula  [21]. 

G/E  composites  have  low  density,  high  stiffness,  and  low  CTE,  making 
them  prime  candidates  for  the  use  in  LSS.  Cammahort,  Rennhack,  and  Coons 
[22]  have  studied  the  effects  of  a  thermal  cycling  enviroment  on  the  dimen¬ 
sional  stability  and  microstructural  integrity  of  G/E  composites.  They 
subjected  their  specimen  to  25  thermal  cycles  between  the  temperature  of 
liquid  nitrogen  (77K)  and  that  of  boiling  water  (373K). 

Bowles,  Post,  Herakovich,  and  Tenney  [23]  and  Bowles  and  Tenney  [24], 
describe  an  experimental  program  conducted  to  determine  the  feasibility  of 
applying  moire  interferometry  to  the  measurement  of  thermally  induced  strains 
in  composites.  They  have  developed  techniques  for  high  sensitivity  moire 
interferometry  by  reflection,  using  a  real  reference  grating  of  1200  lines/mm. 
They  have  determined  the  CTE  for  four  G/E  laminates  in  the  temperature  range 
of  297K  to  422K:  (1)  [0];  (2)  [90];  (3)  [0/+45/90];  and  (4)  [0/90±45]. 

The  dimensional  stability  of  materials  used  in  LSS  must  be  established 
at  both  low  and  high  temperatures.  The  moire  interferometry  system  of  ref¬ 
erences  23  and  24  is  limited  to  making  CTE  measurements  above  room -tempera¬ 
ture.  Short,  Hyer,  Bowles,  and  Tompkins  [25]  have  developed  an  interfer¬ 
ometer  technique  to  measure  small  thermal  strains  associated  with  G/E  com¬ 
posites  over  the  temperature  range  of  116K  to  366K. 

When  working  with  composite  materials,  the  designer  must  be  concerned 
with  the  effects  of  moisture  absorption/desorption.  Shen  and  Springer  [26] 
have  developed  expressions  for  the  moisture  content  of  composite  materials 
as  a  function  of  time.  The  results  apply  when  both  moisture  content  and 
the  temperature  of  the  environment  are  constant.  Bergman  and  Dill  [27J  have 


8 


fctJLii  JLm  ■Kwh  Jh> 


%**  V 


s 


.-v 


n 


\ 

k 

u 


*5. 


V 


i 


■: 

*. 


•3 


conducted  experiments  to  determine  the  effect  of  absorbed  moisture  on  the 
strength  and  stiffness  properties  of  G/E  composites.  Zigrang  and  Bergman 

[28]  have  studied  the  response  of  G/E  sandwich  panels  to  moisture  and 
temperature  transients.  The  application  prompting  this  study  was  the  use 
of  G/E  honeycomb  sandwich  construction  for  the  Shuttle  payload  bay  doors. 

Durability  of  composites  in  the  GEO  environment  is  considered  to  be  the 
principle  uncertainty  associated  with  the  long  term  [25-30  years)  use  of 
composites  in  thin-gage  minimum  weight  LSS.  Tenney,  Slemp,  Long,  and  Sykes 

[29]  describe  two  programs  conducted  at  LaRC  in  FY-79  to  assess  the  radiation 
stability  of  current  composites.  The  goal  of  one  program  is  to  determine  the 
radiation  damage  mechanisms  of  resin  matrix  composites  and  formulate  new 
polymer  matrices  that  are  inherently  more  stable  in  the  space  environment. 

The  thrust  of  the  second  program  is  to  develop  the  materials  technology 
required  for  confident  design  of  LSS  such  as  antennas  and  platforms.  Another 
long  term  concern  for  materials  exposed  to  conditions  found  in  space  is  stress 
corrosion  cracking.  Stress  corrosion  cracking  is  the  result  of  both  a 
sustained  tensile  stress  and  a  corrosive  environment  such  that  a  material 
fails  prematurely.  Marshall  Space  Flight  Center  [30]  has  investigated  this 
problem  and  offers  criteria  for  material  design. 

Because  most  polymeric  materials  are  known  to  undergo  changes -in  physical 
and/or  mechanical  properties  when  exposed  to  electron  and  proton  radiation, 
a  program  designed  to  evaluate  the  effects  of  this  radiation  on  selected 
candidate  composite  materials  has  been  conducted  at  LaRC.  Slemp  and  Santos 
[31]  present  some  of  the  initial  significant  results  from  this  research  program. 
Fornes,  Memory,  Gilbert,  and  Long  [32]  present  data  gathered  during  experimenta¬ 
tion  conducted  to  evaluate  the  effects  of  radiation  on  the  mechanical  properties 
of  epoxy  based  structural  materials  and  to  measure  the  fundamental  radiation- 


generated  events  which  cause  the  changes.  G/E  composite  specimens  and  epoxy 
resins  were  exposed  to  electron  and  gamma  radiation,  followed  by  mechanical 
property  and  fundamental  measurements.  Radiation  induced  changes  in  a 
structural  material  are  initiated  by  energetic  charged  particle  impact  which 
transfers  kinetic  energy  of  the  particle  via  its  coulomb  field  into  electronic 
excitation  energy  of  the  material.  Kamorator,  Wilson,  Chang,  and  Xu  [33]  have 
investigated  the  energy  absorption  and  events  immediately  following  impact  in 
an  effort  to  arrive  at  an  understanding  of  the  radiation  degradation  of 
polymers . 

For  organic  matrix  composites,  microcracking  is  one  of  the  primary  fac¬ 
tors  controlling  the  dimensional  stability  of  the  material.  Permanent  changes 
in  the  thermoelastic  properties  and/or  permanent  residual  strains  can  result 
from  this  type  of  damage.  Bowles  [34]  has  studied  the  effect  of  microcracking 
on  the  CTE  of  G/E  composites. 

Hoggatt  and  Kushner  [35]  discuss  a  thermoplastic  resin  (polysulfane) 
and  compare  the  structural  and  environmental  properties  and  the  fabrication 
and  repairability  of  the  thermoplastic  composite  with  a  typical  epoxy  com¬ 
posite.  It  is  their  conclusion  that  low  labor  costs  exhibited  by  the  thermo¬ 
plastic  composites  make  them  a  priority  consideration  for  use  in  space  struc¬ 
tures.  Fager  [36]  gives  a  good  argument  supporting  the  contention  that  G/E 
is  the  material  to  use  in  LSS. 

A  class  of  materials  which  have  only  recently  become  practical  to  use 
in  structures  are  metal-matrix  composites.  Some  advantages  of  these  materials 
are  better  response  at  high  temperatures  and  the  possibility  of  using  a  matrix 
with  high  ductility  to  increase  fracture  toughness.  Vinson,  and  Chow  [37], 
and  Renton  [38],  as  well  as  Herakowich,  Davis,  and  Dexter  [39]  discuss  various 
applications  of  metal-matrix  composites  along  with  some  material  properties. 


Divecha,  Fishman,  and  Karmarkar  [40]  present  production  techniques, 
material  properties,  and  practical  application  of  Silicon  Carbide/Aluminum 
Composites.  It  is  necessary  to  develop  a  constitutive  damage  model  of 
material  behavior  if  a  structure  is  to  be  designed  to  function  under 
degraded  material  conditions.  Much  work  has  been  done  involving  fracture 
mechanics  in  composite  materials  in  general  [41-48].  Yeow,  et  al.  [49-51] 
have  investigated  the  fracture  behavior  of  Graphite/Epoxy  composites.  Hahn 
[52-55]  has  investigated  the  fracture  rsponse  of  both  polymeric  and  metal- 
matrix  composites. 

Some  of  the  earliest  work  on  modeling  materials  was  for  metals.  Bodner 
et  al.  [56-59],  and  Krieg,  Swearengen,  and  Rohde  [60]  have  constructed 
microphenomenologically-based  models  for  the  response  of  metals. 

Several  damage  models  are  currently  available  for  use  in  composite  materials. 
Maryolin  [61]  has  done  work  on  a  constitutive  model  based  on  a  prescribed 
distribution  of  penny  shaped  cracks.  A  damage  field  is  constructed  by  the 
use  of  the  crack  distribution  and  crack  radius.  With  this  model,  the  effective 
moduli  of  the  material  can  be  calculated.  Laws,  Dvorak,  and  Hejazi  [62]  have 
modelled  a  unidirectional  continuous  fiber  composite  containing  longitudinal 
slit  cracks.  By  assuming  different  magnitudes  of  cracks,  both  a  two-phase 
and  three-phase  model  results.  Horii  and  Nemat-Nasser  [63]  have  developed 
a  model  involving  randomly  distributed  cracks  and  frictional  sliding.  Both 
Margolin,  Dvorak,  and  Nemat-Nasser  use  the  self-consistent  scheme  (i.e., 

Eshelby's  solution)  in  their  models  for  the  determination  of  the  stiffness 
reduction. 

A  plane-stress  continuum  model  for  creep  analysis  of  a  unidirectional 
metal-matrix  composite  was  developed  by  Min  and  Crossman  [64].  This  model 
used  separate  transverse  and  shear  damage  parameters  to  account  for  damage 


caused  by  inhomogeneous  local  distribution  of  stress  in  the  material. 

Hashin  and  others  [65-66]  have  modelled  damage  by  the  use  of  damage 
surfaces.  These  surfaces  are  based  on  S-N  curves  where  N  is  the  number  of 
cycles  to  give  fracture  at  a  given  stress  S. 

Chou  and  Croman  [67]  formulated  a  model  using  the  equal  rank  assumption 
first  introduced  by  Hahn  and  Kim  [68].  This  model  assumes  that  fatigue  life 
distribution  has  the  same  rank  as  in  the  static  case.  The  assumption  was 
based  on  an  overall  decrease  in  residual  strength,  rather  than  specific  cracks. 

Fonseka,  and  Krajcinovic  [69,70]  have  developed  a  thermodynamically-based 
damage  model  for  use  on  brittle  materials  with  continuous  damage.  This  model 
uses  internal  state  variables.  Another  model  using  internal  state  variables 
was  developed  by  Talreja  [71-72].  Talreja's  model  is  a  continuum  mechanics 
damage  model  which  has  damage  from  any  loading  condition  expressed  in  terms  of 
vector  damage  parameters.  Allen  [73]  has  proposed  a  damage  model  for  composite 
materials  based  on  a  continuum/thermodynamics  approach.  In  his  model  both 
local  damage  as  well  as  global  constitutive  equations  are  averaged  over  the 
volume  of  the  laminate.  The  internal  state  variables  of  damage  that  result  are 
vector  valued  history  dependent  functions. 

Schapery  [74-76]  has  constructed  a  model  for  damage  in  a  composite  based 
on  fracture  mechanics  and  nonlinear  viscoelasticity.  This  model  has  been 
primarily  formulated  for  isotropic  particulate  composite  materials,  but  can 
be  extended  to  model  continuous  fiber  composites. 

Research  has  been  performed  at  General  Dynamics  (Forth  Work)  [77]  with 
damage  modeling.  Work  has  been  done  on  a  damage  program  outline  and  a  pre¬ 
liminary  damage  model  exists  at  this  time. 


IV.  Solution  Techniques  for  LSS 

Hie  solution  techniques  for  LSS  can  be  classified  into  several  categories: 
Noor  [78]  classifies  them  into  four  categories:  (1)  Finite  element  modelling, 

(2)  continuum  lattice  modelling,  (3)  discrete  field  method,  and  (4)  periodic 
structure  approach.  This  classification  is  based  upon  repetitive  lattice 
trusses,  which  are  expected  to  be  the  primary  candidates  for  most  of  the  large 
space  structures.  In  an  earlier  status  report  on  latticed  structures  [79] 
the  analysis  methods  are  divided  into  six  categories:  (1)  discrete  field 
analysis,  (2)  space  frame  analysis,  (3)  equivalent  continuum  analysis,  (4) 
nonlinear  analysis,  (5)  experimental  analysis,  and  (6)  dynamic  analysis. 

The  paper  by  Roussos,  Hyer  and  Thornton  [80]  develops  a  procedure  by 
which  many  different  damping  functions  (linear  or  nonlinear)  may  be  incorporated 
in  a  structure.  In  the  finite  element  formulation,  a  solution  technique  is 
developed  by  modifying  the  Newmark  method  for  solving  the  set  of  nonlinear 
equations.  Mehaney,  Thornton  and  Dechaumphai  [81,82]  have  developed  integrated 
Thermal-Structural  finite  element  analysis.  In  this  integrated  approach, 
thermal  and  structural  finite  elements  are  formulated  with  a  common  geometric 
discretization  and  the  thermal  and  structural  finite  elements  are  fully 
compatible. 

Adelman  and  Shore  [83]  describe  the  thermal  analysis  considerations 
for  large  space  structures.  The  paper  by  A.  K.  Noor  [84]  summarizes  recent 
developments  in  the  application  of  reduction  methods  to  nonlinear  structural 
mechanics  problems.  Some  of  the  aspects  of  reduction  methods  discussed  include 
selection  of  basis  vectors  in  nonlinear  static  and  dynamic  problems  and  the 
use  of  reduction  methods  in  conjunction  with  mixed  finite  element  models. 

A  similar  paper  by  Noor  and  Peters  [85]  describes  a  reduced  basis  technique 

and  computational  algorithm  for  predicting  nonlinear  static  response  of  structures. 
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Soni  and  Bogner  [86]  present  a  computer  program  MAGNA-D  for  predicting 
the  response  of  damped  structures  to  steady  state  inputs.  The  analysis 
includes  the  damping  effects  of  viscoelastic  materials  characterized  by 
complex  moduli  and  of  coulomb  friction  at  sliding  interfaces. 

One  of  the  earlier  papers  on  continuum  modelling  is  by  Heki  and  Saka  [87]. 
This  paper  presents  the  stress  analysis  of  lattice  plates  as  anisotrpic  continuum 
plates.  In  this  approach  effective  rigidity  of  lattice  plates  have  been  derived 
for  bending,  stretching,  shearing  and  coupling  rigidity  between  them.  Flower 
and  Schmidt  [88]  present  an  analysis  of  space  truss  as  equivalent  plate.  In 
this  analysis  it  is  assumed  that  all  truss  members  behave  linearly  classically 
and  are  pin-pointed.  Noor  [89]  presents  a  thermal  analysis  of  double  layered 
grids.  These  double  layers  are  useful  because  of  their  low  cost  light  weight 
and  high  stiffness.  Also  due  to  their  ease  of  packaging,  transporting  and 
assembling  in  space,  double  layered  lattice  trusses  are  considered  as  prime 
candidates  for  large  space  structure  application. 

Nayfeh  and  Hefzy  [90]  present  a  continuum  modeling  technique  of  mechanical 
and  thermal  behavior  of  discrete  large  structures.  This  procedure  can  be  applied 
to  a  wide  variety  of  discrete  structures.  In  another  paper  by  the  same  authors 
[91]  two  models,  cubic  and  octetruss  are  compared  and  their  analysis  indicates 
that  from  the  stiff ness-to-density  ratio  the  cubic  model  is  a  better  candidate 
for  the  large  structures. 

Noor,  Anderson  and  Greene  [92]  have  developed  continuum  models  for  large 
repetitive  beam-  and  platelike  structures  with  arbitrary  configurations  subjected 
to  static,  thermal  and  dynamic  loadings.  This  model  accounts  for  local  effects 
in  the  repeating  elements  (which  were  neglected  in  most  of  the  earlier  models) 
of  the  actual  structure.  Noor  and  Anderson  [93]  present  a  simple  procedure 
for  predicting  the  thermoelastic  and  free  vibration  responses  of  large  repetitive 


beam-like  trusses.  In  this  procedure  the  original  lattice  structure  is  replaced 
by  an  equivalent  continuum  beam  model.  The  equivalent  beam  model  accounts  for 
warping  and  shear  deformation  in  the  plane  of  the  cross-section  and  is 
characterized  by  its  thermoelastic  strain  and  kinetic  energies  from  which  the 
equations  of  motion  and  constitutive  relations  can  be  derived. 

A  simple  procedure  is  presented  by  Noor  and  Weisstein  [94]  for  predicting 
the  buckling  loads  associated  with  general  instability  of  beam-like  lattice 
trusses.  Bazant  and  Christensen  [95]  present  a  micropolar  analogy  to  model 
large  grid  frameworks  under  initial  axial  forces.  This  continuous  approximation 
is  formulated  even  for  the  case  of  variable  member  properties  and  variable 
initial  axial  forces.  Also  the  boundary  disturbances  which  rapidly  decay 
with  distance  from  boundary  is  demonstrated.  The  paper  by  Sun  and  Yang  [96] 
present  a  continuum  approach  toward  dynamics  if  gridworks.  In  this  approach 
equations  of  motion,  boundary  conditions  and  constitutive  equations  are  developed 
which  can  be  used  to  analyze  dynamic  and  static  responses  of  gridwork  consisting 
of  orthogonally  interesting  beams  rigidly  connected  at  the  joints.  The  transverse 
vibration  and  wave  propagation  are  investigated  by  continuum  theory.  Comparison 
with  finite-element  solutions  show  that  results  obtained  are  accurate  especially 
in  the  range  of  long  wavelengths. 

Noor  and  Nemeth  [97]  have  developed  micropolar  models  for  the  static,  free 
vibration  and  buckling  analysis  of  repetitive  spatial  beam-like  lattices  with 
rigid  joints. 

The  discrete  field  method  takes  advantage  of  the  regularity  of  the  lattice 
and  is  based  on  construction  of  the  equilibrium  and  compatibility  equations  at 
a  typical  joint  and  then  using  a  Taylor  Series  expansion  to  replace  these 
equations  by  difference  equations.  Renton  [98-101]  has  analyzed  fluxural- 
torsional  behavior  of  grillages,  triangular  mesh  grillages,  space  grids  and 


on  the  gridwork  analogy  for  plates.  Wah  [102]  outlines  a  procedure  for 
analyzing  laterally  loaded  gridworks. 

The  periodic  structure  approach  is  generally  based  on  either:  (a)  the 
combined  use  of  finite  elements  and  transfer  matrix  methods,  which  is  efficient 
only  for  rotationally  symmetric  or  simple  geometries,  or  (b)  the  exact  representa¬ 
tion  of  the  stiff  of  an  individual  member  from  which  the  analysis  of  beam-like 
lattices  with  simple  supported  edges  cam  be  performed.  McDaniel  and  Chang 
[103]  employ  a  finite  element  transfer  matrix  method  to  model  dynamics  of 
rotationally  periodic  large  space  structures.  In  this  approach,  the  finite 
element  transfer  matrix  method  is  employed  to  eliminate  internal  degrees  of 
freedom  from  the  basic  unit  of  rotationally  periodic  large  space  structures. 
Eigenfunctions  of  the  resulting  periodic  unit  transfer  matrix  are  used  to 
obtain  frequency  responses  of  the  complete  structure  without  increasing  the 
analysis  variables.  Interpolation  procedures  are  developed  which  significantly 
reduce  the  required  computations.  The  dimension  of  the  transfer  matrix,  and 
the  number  of  eigenvalues/eigenvectors  extractions  required  in  a  given 
frequency  range. 

Anderson  [104-106]  has  developed  equations  for  general  lattice  structures 
having  repetitive  geometry  for  vibration  and  buckling  analysis.  This  approach 
is  based  on  represnting  each  member  of  the  structure  with  the  exact  dynamic 
stiffness  matrix  and  taking  advantage  of  the  repetitive  geometry  to  obtain  an 
eigenvalue  problem  involving  the  degrees-of-f reedom  at  a  single  node  in  the 
lattice.  According  to  Noor  the  limitations  of  the  periodic  structure  approach 
can  be  removed  by  combining  this  approach  with  the  continuum  approach. 

V.  Damping  Technology 

Much  research  has  gone  into  the  area  of  constrained  layer  damping.  Carne 
[107]  examines  constrained  layer  damping  by  using  finite  element  analysis. 


In  this  paper  he  gives  a  brief  outline  of  the  early  development  (1959  to  1973) 


of  constrained  layer  damping  technology.  The  first  important  paper  on  constrained 
layer  damping  was  by  Kerwin  [108].  In  it  he  develops  theoretical  expressions 
for  the  damping  effectiveness  of  a  constrained  damping  layer  based  on  the  shear 
energy  -  loss  mechanism.  This  was  followed  by  a  more  extensive  paper  by  Ross, 
Ungar,  and  Kerwin.  Using  their  theory,  the  loss  coefficient  and  the  increase 
in  dynamic  stiffness  due  to  the  constrained  layer  could  be  determined.  Their 
theory  is  applicable  for  beams  of  infinite  length  or  finite  beams  of  infinite 
length  or  finite  beams  with  simply  supported  boundary  conditions.  Ungar  [109], 
using  an  energy  method,  was  able  to  calculate  the  loss  coefficient  of  more 
complicated  beam  structures.  Yin,  Kelly,  and  Barry  [110]  verified  Ross,  Ungar, 
and  Kerwin.  Their  paper  dealt  with  the  constrained-layer  damping  mechanism 
for  plates,  beams,  and  tubular  structures.  DiTaranto  [111]  derived  a  6th 
order  PDE  to  describe  the  motion  of  sandwich  beams,  thus  allowing  boundary 
conditions  other  than  simply  supported  for  finite  length  beams.  DiTaranto 
and  Blasingame  [112]  applied  the  PDE  to  obtain  generalized  results  for  the 
composite  loss  factor  and  natural  frequency  of  a  sandwich  beam.  DiTaranto 
derived  his  6th  order  PDE  in  terras  of  longitudinal  displacements.  Mead  and 
Markus  [113]  rederived  the  equation  based  on  lateral  displacements,  and  then 
used  it  [114]  to  determine  the  characteirstic  equation  for  the  resonant 
frequency,  loss  factor  and  modal  root.  Yan  and  Dowell  [115]  derived  a  simpler 
4th  order  PDE  for  the  vibrating  beam  with  constrained  layer  damping,  and  they 
published  [116]  some  experimental  results  in  verification  of  the  theory.  Mead 
[117]  discusses  some  questions  he  has  concerning  the  4th  order  PDE.  Carne's 
paper  was  written  in  1975.  In  1978  DiTaranto  published  "Summary  of  Constrained 
Damping  Literature"  [118].  This  paper  includes  the  literature  survey  of  B.C. 

Nakra  and  has  well  over  100  papers  referenced.  More  recent  papers  on  constrained 


damping  include  Rao  [ll9]  which  uses  an  energy  method,  Johnson  and  Kienholz 
[120],  and  Bogner  and  Brockman  [111]  use  the  finite  element  method.  Soovere 
[122]  wrote  of  a  program  for  the  development  of  dampin'  applicable  to  stringers 
and  frames  on  a  wide-bodied  transport.  The  use  of  constrained  layer  damping 
has  been  extended  to  rings  and  cylinders  by  DiTaranto  [123]  and  Lu,  Douglas, 
and  Thomas  [124]. 

Damping  materials  and  systems  are  represented  herein  by  three  papers.  Jones 
[125]  wrote  an  introduction  to  damping  materials  and  system.  In  it  he  lists 
32  references.  O’Kiefe  [126]  discussed  application  techniques  and  utilities 
for  several  types  of  damping  systems.  He  described  the  microscopic  mechanics 
of  damping  as  well  as  some  methods  for  improving  performance  over  wider  temperature 
and  frequency  ranges.  Gibson  [127]  investigates  the  internal  damping  and  dynamic 
stiffness  of  the  G/E  composite  column  of  reference  160. 

Experimental  investigations  of  damping  were  conducted  by  Plunkett  [128,129] 
and  by  Plunkett  and  Sax  [130]. 

A  few  papers  are  included  which  are  concerned  with  designing  for  vibration 
control.  Rogers  [131]  wrote  an  outline  guide  for  the  design  of  vibration  damping. 
Paul  [132]  illustrates  damping  thermal  design  methodology  by  discussion  of  two 
recent  damping  design  projects.  These  concern  the  TF30P100  inlet  guide  vane 
and  the  B-l  aft  equipment  bay. 

The  following  papers  are  directed  toward  active  control:  Horner  [133], 
Herzberg,  Johansen,  and  Stroud  [134],  Skelton  [135,136],  and  Montgomery  [137]. 

Snowden  [138]  is  an  often  referenced  text  on  the  subject  of  vibration  and 
damping . 

Foss  [139]  derived  orthogonality  relationships  between  the  eigenvectors 
of  damped  linear  dynamic  systems  with  lumped  parameters,  and  was  able  to  establish 
coordinates  which  uncouple  the  equations  of  motion. 


VI.  Preliminary  Design/Analysis/Experimental  Investigation 

Any  LSS  will  be  built  using  the  STS  and  must  therefore  be  designed  for 
compatibility  with  this  system.  Roebuck  [140]  describes  shuttle  considerations 
for  the  design  of  LSS.  Culbertson  and  Bold  [141]  references  many  Shuttle  user- 
oriented  documents.  Christensen  [142]  requests  input  necessary  to  develop  LSS 
design  guidelines  and  standards. 

There  are  basically  three  methods  of  placing  large  structures  in  space: 

(1)  automated  fabrication,  wherein  a  manufacturing  system  is  placed  in  orbit 
which  when  supplied  with  raw  materials  will  build  truss  like  structures;  (2) 
deployable  structures  built  on  Earth  and  collapsed  into  a  compace  transportation 
configuration  to  be  unfurled  in  space;  and  (3)  erectable  structures  built  from 
pieces  carried  into  orbit.  Regardless  of  the  method  of  construction,  these  will 
be  first  deployed  in  LEO  and  then  transferred  to  GEO. 

Powell  [143],  and  Powell  and  Browning  [144]  describe  a  baseline  LSS  concept 
developed  for  the  NASA-JSC,  SCAFE.  A  beam  builder  and  raw  materials  are  placed 
in  a  566  km  circular  orbit  by  the  STS.  The  beam  builder,  working  from  the  Shuttle 
automatically  fabricates  four  triangular  beams,  each  200m  in  length.  The  beam 
builder  then  moves  along  the  length  of  these  to  build  and  attach  shorter 
cross  beams.  The  structure  is  to  be  built,  and  preliminary  tests  run,  in  a  7 
day  mission.  Slysh  and  Kugath  [145]  describe  a  LSS  automated  assembly  technique 
(LSAT)  developed  to  provide  a  starting  point  for  a  feasibility  examination 
of  fully  mechanized  space  structures  assembly  and  maintenance.  The  LSAT  concept 
is  applicable  to  structures  ranging  in  size  from  hundreds  to  tens  of  thousands 
of  meters.  Muench  [146]  describes  an  aluminum  beam  builder  which  has  been 
built  and  ground  tested  by  Grumman  Aerospace.  It  could  be  used  to  produce 
the  components  necessary  for  assembling  LSS.  Goodwin  [147]  discusses  the  concept 
of  a  space  platform  to  be  used  as  a  manufacturing  base  for  the  beam  builder  of 


reference  100.  The  platform  could  house  solar  arrays,  radiators,  space  construc¬ 
tion  equipment,  laboratories,  housing  modules,  etc. 

The  second  means  of  placing  LSS  in  orbit  is  the  deployable  method.  These 
structures  are  assembled  on  Earth,  packaged  in  a  compact  collapsed  configuration 
to  be  unfurled  in  space.  This  type  structure  can  never  reach  the  dimensions 
necessary  for  many  of  the  proposed  applications  of  LSS,  but  can  be  utilized 
for  larger  aperture  antennas.  Fager  and  Garriott  [148]  describe  a  concept 
capable  of  deploying  a  reflector  up  to  91.5  M  in  diameter.  Spring  loaded  trusses 
retract  for  transportation  and  are  released  in  space.  A  30.5  M  diameter  antenna 
can  be  packaged  into  a  3.05  M  payload  envelope.  Deployment  is  automated. 
Chadwick  and  Woods  [149]  discuss  the  work  being  done  at  LMSCO  to  expand  the 
maximum  size  of  deployable  antennas.  Their  work  is  based  on  the  proved  WRAP- 
RIB  design.  With  new  materials  and  design  improvements,  the  WRAP-RIB  structure 
can  now  be  extended  to  support  antennas  of  over  200  M  in  diameter.  These 
designs  can  be  packaged  for  ascent  on  the  STS  and  are  self  deploying  in  space. 

A  hybrid  method  connecting  the  erectable  concept  with  the  deployable 
concept  proposes  the  construction  of  LSS  by  the  mating  of  deployable  modules. 

Agan  [ 150]  contributes  to  the  overall  technical  data  base  for  various  programs 
using  deployable  structures.  His  paper  is  generic  in  nature,  intended  to 
provide  general  information  and  data  to  emerging  space  programs  employing 
erectable/deployable  structural  concepts.  His  basic  structural  building 
block  is  a  double  cell  double  folding  cubic  model  (3M  x  3M  in  cross-section 
and  6M  long).  Heartquist  [151]  relates  the  results  of  twin  contracts  awarded 
to  General  Dynamics-Convair  and  Martin  Marietta  to  develop  preliminary  design 
of  a  183M  truss-like  spacecraft.  Armstrong,  Skoumal,  and  Straaver  [152]  have 
conducted  a  building  block  structures  study  for  LSS.  Based  on  proposed  uses  of 
LSS,  they  decided  that  a  modular  planar  truss  structure  and  a  long  slender 
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boom  concept  would  serve  as  building  block  approaches.  The  structural  configura¬ 
tions  chosen  for  the  study  were  the  tetrahedral  planar  turss  and  lattice  column. 
Articulating  joints  and  collapsing  details  are  defined  to  allow  packaging 
densities  to  approach  values  that  effectively  utilize  the  payload  mass  capability 
of  the  Shuttle.  A  deployable  planar  truss  building  block  module  is  described 
in  parametric  form  to  allow  a  broad  application  base  for  this  type  of  structure. 
Barclay,  Brogren,  and  Skoumal  [153]  also  deal  with  the  subject  of  reference 
152.  A  scenario  is  developed  in  [153]  which  would  enable  the  placing  of  a 
300M  by  300M  platform  in  LEO  using  the  STS.  The  platform  is  assembled  using  7 
deployable  tetrahedral  modules  [120M  by  104M  in  size)  mounted  to  an  erectable 
truss  structure  consisting  of  102  columns.  The  structural  response  was  determined 
for:  (1)  a  simple  open  truss  subjected  to  the  natural  thermal  radiation,  and  (2) 

a  more  complex  truss  subjected  to  on-board  heating  as  well  as  the  natural  thermal 
radiation. 

Erectable  concepts  deal  with  the  building  of  LSS  from  beam  or  column  elements 
transported  to  LEO  by  the  STS.  Much  work  has  been  done  in  this  area  at  LaRC 
by  Heard,  Bush,  Walz,  Card,  Mikulas,  et  al .  This  work  has  dealt  primarily 
with  the  tetrahedral  truss  configuration.  Mikulas,  Bush,  and  Card  [154] 
have  written  a  good  paper  on  the  tetrahedral  truss.  In  the  paper,  the  truss 
is  defined,  analyzed  (using  a  continuum  method),  and  tested  experimentally. 

The  tetrahedral  truss  is  an  excellent  configuration  for  preliminary  studies 
because  of  its  psuedo-isotropic  elastic  properties  and  because  it  is  con¬ 
structed  of  all  identical  column  members.  Also,  work  done  by  Boeing  identified 
the  tetrahedral  truss  as  having  the  lightest  weight  and  highest  first  mode 
frequency  when  compared  to  other  generic  trusses  of  equal  depth  and  planform 
area.  Reference  154  considers  both  aluminum  and  G/E  tubular  columns.  Work 
by  Bush,  Mikulas,  and  Heard  [155]  is  very  similar  to  reterence  154.  Heard,  Bush, 


Walz,  and  Rehder  [156]  have  conducted  optimization  studies  to  examine  minimum 
mass  structural  proportions  of  deployable  and  erectable  tetrahedral  truss 
platforms  subjected  to  the  integrated  effects  of  various  design  requirements. 
Considerations  integrated  into  the  design  process  are:  (1)  lowest  natural 
frequencies  of  the  platform  and  individual  platform  components;  (2)  packaging 
constraints  imposed  by  the  Shuttle  cargo  bay  capacity;  (3)  initial  curvature 
of  the  struts;  (4)  column  buckling  of  the  struts  due  to  gravity  gradient, 
orbital  transfer,  strut  length  tolerance,  or  design  loads;  and  (5)  lower  limits 
for  strut  diameter  and  wall  thickness.  Katz  and  Pankopf  [157]  describe  a  three 
step  procedure  for  assembling  LSS.  They  use  the  tetrahedral  truss  and  the 
nestable  tapered  column  of  reference  159. 

Heard,  Bush,  and  Walz  [158]  have  tried  to  determine  the  size  of  LSS  that 
is  best  suited  to  a  deployable  concept  versus  an  erectable  concept.  They  used 
the  tetrahedral  truss  for  mathematical  modelling  purposes.  Preliminary  orbital 
transfer  investigations  indicate  that  deployable  platforms  of  up  to  200M  span 
may  be  placed  in  GEO  with  a  single  Shuttle  flight  using  a  constant  thrust 
chemical  propulsion  system  which  limits  initial  accelration  to  O.Olg  or  less. 

Park  and  Winget  [159]  wrote  the  only  paper,  found  by  this  author,  which 
proposed  the  use  of  nonperiodic  LSS.  The  design  and  analysis  of  LSS  for  dynamic 
loads  requires  the  consideration  of  wave  propagation,  transient  response,  and 
steady-state  vibration  problems.  The  desirable  intrinsic  structural  charac¬ 
teristics  of  these  problems  are  good  dispersive  properties,  rapid  decay  of  the 
transients,  and  an  optimum  distribution  of  frequency  spectrum.  Periodic 
trusses  such  as  the  tetrahedral  truss  or  the  octetruss  have  two  major  drawbacks. 
First,  periodic  trusses  can  be  considered  to  be  well  tuned,  and  therefore  they 
do  not  have  good  wave  dispersion  characteristics.  Second,  if  the  dimensions 
of  the  platform  are  fixed  then  the  frequency  spectrum  of  the  structure  is  also 


fixed,  thus  leaving  no  room  for  frequency  modification  other  than  through 
redesign  and/or  effective  vibration  control  devices.  Park  and  Winget  propose 
improving  wave  dispersion  and  controlling  the  frequency  spectrum  by  the  internal 
rearrangements  of  nonperiodic  lattices. 

Preliminary  studies  indicate  that  LSS  will  be  large  area  low  mass  truss 
type  structures  constructed  of  long  slender  lightly  loaded  compression  elements 
may  most  likely  of  G/E.  Bush  and  Mikulas  [160]  have  developed  a  nestable  tapered 
column  to  be  used  as  the  structural  element  of  the  tetrahedral  truss  developed 
at  LaRC.  The  nestable  tapered  column  permits  achievement  of  weight  critical 
as  opposed  to  volume  critical  payloads  for  the  STS.  Structures  using  the 
nestable  tapered  column  will  have  to  be  assembled  in  orbit  eliminating  the 
use  of  deployable  modules.  It  is,  however,  the  opinion  of  the  authors  that 
for  very  large  (km  size)  structures,  the  gain  in  mass  placed  in  orbit  per  launch 
makes  the  nestable  tapered  column  feasible.  The  nestable  tapered  column  concept 
will  require  a  "break  or  center  joint"  in  the  middle  of  each  truss  element  and 
a  cluster  joint  for  each  node.  Heard,  Bush,  and  Agranoff  [161]  have  investigated 
the  effect  of  the  mid-column  break  on  stability  behavior.  Since  the  columns  will 
function  in  an  extremely  low  strain  range  they  have  also  investigated  the  validity 
of  using,  for  analysis  predictions,  standard  material  properties  for  G/E  determined 
from  coupons  tested  udner  much  larger  strain  ranges.  Some  preliminary  buckling 
test  results  for  the  G/E  nestable  tapered  column  appear  in  reference  109.  The 
purpose  of  reference  161  is  to  present  detailed  data  on  the  columns  and  center 
joint  of  the  nestable  tapered  columns  concept,  and  to  present  buckling  data 
for  additional  columns  as  well  as  a  tripod  arrangement  of  these  columns  using 
a  cluster  joint.  Gibson  [127]  has  investigated  the  damping  characteristics  and 
the  dynamic  stiffness  of  the  G/E  composite  materials  used  in  the  restable  tapered 
column.  He  has  also  compared  the  relationships  between  dynamic  properties 


derived  from  tests  of  small  specimens  and  those  of  full  scale  columns.  It  is 
his  conclusion  that  small  specimens  are  representative.  Mikulas  [162]  is 
concerned  primarily  with  design  methods  for  long  lightly  loaded  columns,  and 
with  a  technique  for  comparing  the  masses  of  different  column  concepts.  He 
looks  at  four  column  concepts:  (1)  a  tubular  column,  (2)  a  three  longeron 
truss  column  constructed  from  tubular  members,  (3)  a  three  longeron  truss 
column  constructed  from  solid  rod  members,  and  (4)  a  tubular  column  with  open 
gridwork  "isogrid"  walls.  Optimization  of  LSS  is  generally  considered  in 
terms  of  minimum  weight.  Yoo  [163]  considers  the  optimization  of  a  three 
longeron  truss  column  constructed  from  tubular  members  (one  of  the  columns 
that  Mikulas  looked  at).  The  optimization  procedures  are  based  on  designing 
for  a  column  with  initial  imperfections.  Yoo  compares  his  results  to  those 
of  Mikulas.  Many  LSS  make  use  of  the  low  mass  and  deployability  of  pretensioned 
structures  to  achieve  efficient  designs.  These  structures  use  tension  elements 
(cables,  rods,  and  membranes)  to  provide  stiffness  and  stability  of  structural 
systems.  To  understand  the  fundamental  structural  characteristics  of  pretensioned 
structures,  Belvin  [164]  has  performed  analyses  and  tests  of  some  simple  configu¬ 
rations.  The  buckling  and  vibration  behavior  of  a  pretensioned  stayed  column 
are  discussed  in  detail.  Crawford  and  Benton  [163]  have  determined  the  axial 
compressive  strengths  for  lattice  columns  that  are  initially  wavy  in  both  local 
and  overall  modes.  The  compressive  strength  is  found  to  be  less  than  either 
the  classical  local  or  overall  buckling  strengths.  These  strength  reductions 
are  anticipated  to  be  significant  for  lightweight  space  columns  or  large 
slenderness  ratio.  The  study  was  prompted  by  a  need  to  predict  the  axial 
compressive  strength  of  lattice  columns  which  were  being  considered  for  the 
diagonal  spars  of  an  800M  square  solar  sailing  spacecraft. 
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The  largest  application  of  LSS  technology  being  considered  is  the  SPS. 
Studies  have  been  conducted  regarding  satellites  that  would  produce  a  useful 
output  of  from  5  to  10  GW  at  Earth  receiving  stations.  The  major  dimensions 
of  such  a  satellite  will  be  several  kilometers.  For  instance,  Nansen  and 
DiRamio  [166]  describe  the  structure  for  a  SPS  with  the  following  dimensions 
for  the  solar  collector:  21.28  by  5.30  by  0.47  km.  It  employs  660M  by  660M 
square  bays  470M  deep  as  repeating  elements.  Most  of  the  studies  conducted 
thus  far  involve  a  planar  solar  array  with  either  one  or  two  attached  MPTS's. 
The  large  solar  collector  is  a  source  of  electrical  power  (photovoltaic  cells 
are  utilized)  which  is  converted  to  microwave  electromagnetic  energy  and  re¬ 
layed  from  the  collection  site  in  GEO  to  the  ground  receiving  system.  The 
feasibility  of  the  relay  link  depends  on  the  capability  of  keeping  the  beam 
center  of  each  MPTS  pointed  at  the  ground  receiving  antenna  (rectenna),  and 
the  phase  error  over  the  aperture  small  enough  so  that  efficiency  is  high. 
"Achievable  Flatness  in  a  Large  Microwave  Power  Antenna  Study"  [167]  is  the 
structure  of  the  MPTS,  and  (b)  their  contribution  to  slope  error  of  the  antenna 
array.  "Satellite  Power  System,  Concept  Development  and  Evaluation  Program  - 
Reference  System  Report"  [168]  describes  a  reference  system  concept  developed 
as  the  result  of  twin  system  definition  studies  conducted  by  Boeing  Aerospace 
Company  and  Rockwell  International.  The  reference  system  suggests  implemeting 

two  5  GW  SPS  systems  per  year  for  30  years  starting  in  the  year  2000.  The 
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concept  uses  a  planar  solar  array  (about  55  km  )  built  on  a  graphite  fiber 
reinforced  thermoplastic  structure.  A  one  km  diameter  phased  array  MPTS  is 
mounted  on  one  end  of  the  satellite. 

As  the  need  arises  for  larger  antennas  operating  at  higher  frequencies, 
structural  stiffness  becomes  a  key  design  consideration.  Structural  stiffness 
will  significantly  influence  thermal  distortion,  point  accuracy,  and  natural 


frequency.  Fager  [169]  indicates  that  four  design  criteria  must  be  considered 
regarding  stiffness:  contour  control,  thermal  distortion,  ground  test  capability, 
and  most  critical,  pointing  accuracy.  He  discusses  this  general  problem  area 
with  respect  to  the  geodetic  truss.  Analytical  results  indicate  that  a  careful 
selection  of  materials  and  truss  design,  combined  with  accurate  manufacturing 
techniques,  can  result  in  very  accurate  surfaces  for  large  space  antennas. 

The  purpose  of  the  paper  by  Hedgepth  [170]  is  to  investigate  the  influence  of 
the  accuracy  requirements  on  the  design  geometry  of  the  structure  for  four 
general  types  of  structural  configurations:  tetrahedral  truss,  geodesic  dome, 
radial  rib,  and  pretensioned  truss.  Requirements  on  the  panel  sizes  for  a 
mesh  reflection  surface  are  obtained.  The  influence  of  fabrication  imperfections 
on  the  surface  error  is  also  evaluated.  An  initial  evaluation  is  taken  of  the 
effect  of  thermal  stress  on  the  accuracy.  Finally,  an  evaluation  is  made  of 
the  diameter-wavelength  ratios  potentially  feasible  as  limited  by  fabrication 
imperfections  and  thermal  strains.  As  mentiond  already,  one  of  the  sources 
of  error  in  rf  antennas  is  the  lack  of  dimensional  precision  of  the  surface. 
Hedgepeth  [171]  presents  an  approach  for  estimating  the  amount  of  error  caused 
by  random  dimensional  imperfections  of  the  many  structural  elements  which  make 
up  a  truss-type  antenna.  A  principle  of  equivalence  between  the  analyses  of 
statistical  error  and  of  the  natural  vibration  frequencies  of  the  structure  is 
developed . 

Bellagamba  and  Yang  [172]  have  developed  a  technique  for  constrained  parameter 
optimization  and  have  applied  it  to  the  minimum-mass  design  of  truss  structures. 
The  procedure  employs  and  exterior  penalty  function  to  transform  the  constrained 
objective  function  into  an  unconstrained  index  of  performance  which  is  minimized 
by  the  Gauss  method. 

Soon,  very  large  spacecraft  will  be  assembled  in  LEO  and  then  transferred 
to  GEO.  These  lightweight  and  very  flexible  spacecraft  will  require  specially 


designed  propulsive  systems  for  the  orbital  transfer.  Interactions  between 
such  structures  and  the  propulsive  thrust  forces  can  have  significant  influence 
on  the  overall  spacecraft  configuration,  structural  design,  and  control 
functions  to  a  greater  degree  than  experienced  with  more  conventional  spacecraft. 
Kunz  [173]  examines  some  of  these  structural/propulsive  interactions,  deduces 
significant  propulsion  requirements,  and  illustrates  some  concepts  with  design 
characteristics . 

Oren  and  Cox  [174]  provide  an  evaluation  of  heat  rejection  techniques 
applicable  to  multihundred-kW  space  platforms. 

LSS  preliminary  and  conceptual  design  requires  extensive  iterative  analysis 
because  of  the  effects  of  structural,  thermal,  and  control  intercoupling.  LaRC 
is  developing  a  CAD  program  that  will  permit  integrating  and  interfacing  of 
required  LSS  analyses.  The  primary  objective  of  this  program  is  the  implementation 
of  modeling  techniques  and  analysis  algorithms  that  permit  interactive  design 
and  trade-off  studies  of  LSS  concepts.  Farrell  [175]  presents  an  overview 
of  the  status  of  the  program  and  the  capabilities  added  by  Martin  Marietta. 
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Effect  of  Deqradation  of  Material  Properties  on  the 
Dynamic  Response  of  Large  Space  Structures 
by 

S.  Kalyanasundaram,  J.  D.  Lutz,  W.  E.  Haisler,  and  D.  H.  Allen 

ABSTPACT 

In  this  paper  the  effect  of  degradation  of  material  properties  on 
structural  frequencies  and  mode  shapes  of  Laroe  Soace  Structures  is  invest¬ 
igated.  The  difficulty  and  cost  of  maintenance  of  LSS  make  it  a  necessity  to 
desiqn  these  structures  to  operate  with  a  certain  damage.  The  techniques  for 
achieving  modal  control  will  require  a  more  accurate  knowledge  of  modal  char¬ 
acteristics. 

The  sensitivity  studies  conducted  on  selected  space  truss  and  plane  frame 
problems  indicate  that  the  deqradation  of  material  properties  may  have  a 
significant  effect  on  the  structural  mode  shapes  and  frequencies.  For  even 
small  amounts  of  dearadation(lOT) ,  frequencies  and  location  of  nodes  may  change 
siqnificantly.  It  is  clear  that  these'  effects  must  be  taken  into  consideration 
when  designing  the  control  systems  for  large  space  structures. 


INTRODUCTION 


It  is  projected  that  advanced  aerospace  materials  will  be  required  in  the 
future  generation  of  space  structures-  To  make  the  optimum  use  of  these  materials, 
the  effect  of  environmental  conditions  and  complicated  loading  histories  need  to  be 
studied.  One  of  the  important  consequences  of  the  extreme  loading  conditions  will 
be  the  degradation  of  the  n  iterial  properties  of  the  structural  components. 

The  observed  reduction  in  stiffness  is  often  a  good  indication  of  damage  deve¬ 
lopment  in  the  life  of  these  structural  components.  Degradation  of  material  proper¬ 
ties  usually  leads  to  the  reduction  of  stiffness  which  affects  the  dynamic  response. 

In  this  report,  sensitivity  studies  will  be  presented  which  investigate  the  effect  of 
stiffness  on  structural  frequencies  and  mode  shapes. 

IMPORTANCE  OF  STRUCTURAL  FREQUENCIES  AND  MODE  SHAPES 

The  advent  of  the  Space  Shuttle  has  made  possible  the  development  of  Large 
Space  Structures  (LSS).  Control  systems  for  stabilizing  and  maneuvering  these  very 
large  space  structures,  especially  those  for  precise  pointing,  will  require  extensions 
of  our  present  technology. 

Large  size  by  itself  does  not  arouse  concern,  but  the  structural  flexibility 
resulting  from  minimizing  tp,e  structural  weight,  while  increasing  physical  dimen¬ 
sion  may  present  problems.  At  some  point,  deformations  from  structural  flexibility 
result  in  sufficiently  large  amplitudes  and  very  low  frequencies  (.01  to  10Hz) 
which  create  new  complications  for  control  designers. 

As  an  example  of  the  precision  being  sought  [1],  a  typical  radiometry  application 
may  utilize  a  200-m  antenna  with  a  effective  beam  width  of  0.01  deg.  and  have 
requirements  limiting  the  vib*ratory  beam  shift  to  less  than  0.005  deg.  and  dynamic 
surface  distortions  to  less  than  1mm.  Maneuvering  or  maintaining  altitude  of  such 


2 


•  i  satellite  leads  lu  flexible  body  mol  ion  and  ibis  must  be  well  predicted  and  cen¬ 
tre  1  I  ed . 

The  importance  of  interaction  between  the  control  systems  and  vibratory  res¬ 
ponse  has  influenced  alot  of  research  in  LSS  control  system  [2-4].  The  current 
practice  of  guaranteeing  a  large  separation  between  modal  frequencies  and  the  band¬ 
width  of  control  will  not  be  adequate  in  future  applications.  The  combination  of 
large  size  and  payload-weight  restrictions  will  drive  structural  frequencies  down 
and  need  for  more  accurate  pointing  will  drive  the  control  system  bandwidth  up. 

When  the  frequency  separation  becomes  impossible,  there  exists  a  need  for  adaptive 
control  systems.  This  leads  to  further  research  in  design  of  structural  control 
systems  ac tuator sensor  placement  and  distributed  sensing  and  actuator  as  opposed  to 
co-located  sensors  and  actuator  . 

The  techniques  for  achieving  modal  control  will  require  a  more  accurate  know¬ 
ledge  of  modal  characterisctics .  Optimum  sensor  and  actuator  placement  will  be 
greatly  influenced  by  modal  effects  which  must  be  known  to  a  greater  degree  of  pre¬ 
cision. 

The  difficulty  and  cost  of  maintenance  of  LSS  make  it  a  necessity  to  design 
these  structures  to  operate  with  a  certain  degree  of  damage.  Current  research  on 
damage  characterization  of  advanced  aerospace  materials  [5-9]  indicates  a  stiffness 
reduction  as  high  as  25%  that  can  occur  during  the  life  span  of  these  materials.  Thus, 
there  is  clearly  a  need  to  study  the  effect  of  stiffness  reduction  in  structural 
frequencies  and  mode  shapes. 

MODE!  DESCRIPTION 

The  effect  of  stiffness  reduction  (resulting  from  material  damage)  is  investigated 
In  this  report  for  selected  space  truss  and  plane  frame  problems.  The  frequencies  and 
mode  shapes  were  *'ound  using  finite  element  program.-,  incorporating  standard  space 
truss  and  plane  frame  elements. 
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The  damage  is  modeled  in  this  study  as  a  reduction  in  stiffness  for  ail  the  .struc¬ 
tural  elements.  ihe  amount  of  reduction  in  stiffness  is  assumed  to  be  proportional 
to  the  stress  in  the  structural  element  with  a  maximum  reduction  occuring  in  the  high¬ 
est  stressed  element. 

The  following  equations  govern  the  severity  and  distribution  of  damage: 

L*  -E  [  1-kJ  ] 

°  — y 

L*  =  Reduced  modulus  of  the  element 

E  =  undecraded  modulus  of  the  element 

o  ° 

k=  Constant  of  degradation 


3=  stress  in  the  element 


i 


a 


% 

£ 


t; 


Y  =  Yield  stress 


The  linear  form  of  the  material  damage  law  is  used  herein  for  simplicity.  In 
reality,  the  damage  law  will  be  more  complex  and  history  dependent.  Damage  laws  are 
currently  being  developed  [9-11]  and  will  be  considered  in  future  research.  The  effect 
of  residual  stress  is  not  considered  here  and  a  quasi-elastic  behavior  is  assumed. 

In  addition  to  the  above  model,  rand.oui  reduction  in  stiffness  was  also  consi¬ 
dered  to  study  the  effect  of  degradation  of  some  members  that  might  experience  Severe 
loading  histories. 

Figure  [1]  illustrates  the  geometry  of  the  example  space  truss.  This  structure 
has  10  bays  and  is  fixed  at  one  end.  124  space  truss  elements  are  used  in  the  finite 
element  model.  In  the  initial  undegraded  conf  igurat  ion,  the  material  properties  are 
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.same  for  all  members  with  the  following  values: 

Material  type  -  Graphite  Epoxy  (Hexel) 

Young's  modulys  E  =  21.5  x  10^  psi 

2 

Cross  sectional  area  A  =  1.0  in 

3 

Density  p  =  0.065  lb/in 

Figure  [2]  is  the  finite  element  model  of  the  plane  frame.  21  frame  elements 
are  used.  The  undegraded  material  properties  (same  for  all  members)  are  given  below: 

Material  type  -  graphite  epoxy  (Hexel) 

6 

Young's  modulus  E  =  16  x  10  psi 

Cross  sectional  area  A  =  1.37  in^ 

Moment  of  Inertia  I  =  0.033  in’4 

z 

3 

density  p  =  0.06  lb /in 

Figure  [3]  is  the  plane  truss  model  with  the  same  material  properties  as  the 
space  truss  model..  12  Truss  elements  are  used  to  model  this  case  to  study  the  effect 
of  random  reduction  in  stiffness. 

DISCUSSION  OF  RESULTS 

In  figure  [14],  the  first  two  natural  frequencies  are  plotted  for  different  mag¬ 
nitude  of  damage  for  the  space  truss.  The  properties  were  degraded  by  using  stress 
obtained  from  displacements  orresponding  to  the  first  mode  shapes.  Consequently  the 
damage  or  stiffness  degradation  varies  spatially  from  element  to  element  and,  in 
general,  is  greatest  at  the  fixed  end  (:■:  =  ( i  ) 

The  effects  of  damage  on  frequencies  of  the  truss  are  clear.  Increase  of  damage 


results  in  a  decrease  ol  frequencies  due  to  the  reduction  in  the  stiffness. 


A  maximum 


damage  of  25%  results  in  a  frequency  reduction  of  nearLy  10%  for  the  first  two 
natural  frequencies. 

The  vertical  deflection  of  the  free  end  of  the  truss  corresponding  to  the  first 
mode  is  plotted  in  Figure  (5)  for  the  case  of  undegraded  and  25%  degraded  cases. 
Damage  has  little  effect  on  this  mode  shape  with  the  difference  in  deflection  at  any 
node  lit ing  less  than  2%.  The  second  mode  shape  gives  interesting  results.  Figure  (6) 
shows  deflections  for  different  damage.  The  plot  indicates  a  change  in  modal  deflec¬ 
tion  of  9%  for  5%  damage  to  nearly  73%  change  for  25%  damage.  The  location  of  the 
mode  (zero  displacement)  changes  drastically.  Such  a  movement  of  the  systems  mode 
would  require  very  "robust"  control.  The  third  mode  shape  produces  similar  results 
with  large  changes  in  the  modal  displacement  and  is  plotted  in  Figure  7.  Thus  the 
damage  has  a  significant  effect  on  natural  frequencies  and  mode  shapes. 

For  the  study  of  random  damage  effects  on  the  12  element  plane  truss  (Figure  3), 
the  following  three  cases  are  used:  50%  stiffness  reduction  to  element  9  (2)  50% 
stiffness  reduction  to  elements  8,  9,  and  10,  and  (3)  50%  stiffness  reduction  to 
elements  4,  5,  6,  8,  9,  10  and  12.  Figure  (8)  gives  the  plot  of  the  first  two  fre¬ 
quencies  for  different  cases.  The  frequencies  decrease  as  tthe  properties  of  more 
elements  are  degraded.  This  is  due  to  the  reduction  in  stiffness.  Figures  9  and  10 
compare  the  first  two  mode  shapes  for  the  different  cases.  There  are  some  variations 
in  the  mode  shapes  for  the  different  cases. 

The  effect  of  damage  on  the  plane  frame  problem  is  investigated  by  varying  the 
maximum  reduction  of  stiffness  from  5%  to  25%.  It  is  found  that  the  vertical  mem¬ 
bers  and  the  members  near  the  fixed  end  are  more  affected  when  the  properties  were 
degraded  in  a  similar  fashion  to  space  truss  problem  (i.e.,  using  stress  calculated 
from  displacement  corresponding  to  the  first  mode).  Figure  11  illustrates  the  effect 
of  damage  in  the  first  three  frequencies.  The  frequencies  decrease  with  the  increase 
of  damage  from  5%  to  25%.  Figures  12,  13  and  14  give  plots  of  vertical  displacements 
corresponding  to  the  first  tlireee  mode  shapes  along  the  elements  i,  18  and  15.  There 


is  a  small  variation  corresponding  to  the  first  and  tliinl  inode  shapes.  For  the 
second  there  is  significant  variation  in  the  mode  shape  for  damaged  and  undamaged 
cases . 

For  the  study  of  random  damage  effects  on  the  plane  frame,  the  following  three 
cases  are  used:  (l)  25%  stiffness  reduction  to  element  15,  (2)  25%  stiffness  reduc¬ 
tion  to  elements  8,  15  and  (3)  25%  stiffness  reduction  to  elements  8,  15  and  1. 

These  elements  are  located  at  the  top  of  the  frame  and  might  be  damaged  in  this 
fashion  due  to  localized  impactm  heating,  etc.  The  frequencies  decrease  as  stiffness 
of  more  number  of  elements  are  reduced  and  che  first  three  mode  shapes  are  plotted 
in  figures  15,16  and  17  . 

To  investigate  the  effect  of  damage  on  the  transient  response,  the  space  truss 
model  (figure  1.)  is  subjected  to  an  impact  load  of  10,000  lbs  at  nodes  29  and  30 
in  the  transverse  direction.  The  duration  of  the  load  is  one  second.  Figure  18  gives 
the  plot  of  vertical  displacement  vs.  time  for  the  node  41  at  the  free  end.  There  is  a 
change  in  both  the  frequency  and  the  amplitude  of  response. 

The  next  loading  is  a  cycle  temperature  loading  with  a  period  of  one  day.  The 
temperature  distribution  is  uniform  over  the  structure.  Figure  19  depicts  the  ver¬ 
tical  response  of  node  41  vs.  time  for  a  maximum  temperature  difference  of  250°  at 
the  mid  point  of  each  element.  The  amplitude  increases  with  the  loss  of  stiffness. 

The  amplitude  increases  by  17%  for  the  most  severly  damaged  case  (a  maximum  of  25  : 
stiffness  reduction). 

CONCLUSION 

This  study  confirms  the  need  for  a  more  accurate  knowledge  of  modal  characteristics 
for  large  space  structures.  Degradation  of  material  properties  may  have  a  signifi¬ 
cant  effect  on  the  structural  mode  shapes  and  frequencies.  For  even  small  amounts  of 
degradation  (10%),  frequencies  and  location  of  nodes  may  change  significantly.  It  is 
clear  that  these  effects  must  be  t.ken  in t o  CO  ns idera t i on  when  designing  the  control 
systems  for  large  spare  structures. 

in 


FIGURE  NO.  4 


FIRST  MODE  SHAPE 


23.0 
22.0 
21.0 
20.0 
19.0 
18.0 
17.0 
18.0 
15.0 
o  14.0 
^  13.0 

sj  12.0 

1  n.o 

10.0 


■  V  v  V  .■ 


FIGURE  NO.  5 


^  *  .  "  •  •  l  *  -  •  -  -  •  "  -  a-  •  »  •  -  ->  *>  •>  •  *  *  V 


12 


SECOND  MODE  SHAPE 


FIGURE  NO.  6 


DISPL.(/100) 


THIRD  MODE  SHAPE 

rmnoATivn  dvdphvt 

DEGRADED  25  PERCENT 

288. 


HORIZONTAL  POSITION 


FIGURE  NO.  7 


frequency  vs.  damage 


FIGURE  no.  8 


I  I  I  I  I  I  I  I  T  I  1  !!  I  I  I  I  I  I  I  I  PI  T  I  I  I  I  1 

FIRST  MODE  SHAPE 


-  UNDAMAGED 

-  1  B  LB  ME  NT  DAMAGED 

-  3  ELEMENTS  DAMAGED 

-  7  ELEMENTS  DAMAGED 


SECOND  MODE  SHAPE 


FIGURE  NO.  10 


FREQUENCY  (HERTZ) 


INI 


FIGURF  NO.  11 


18 


SECOND  MODE 


. .  ■  UNDAMAGED 

— — — -  DBGRADHD  5  PURGE NT 

- - -  DBGRADHD  15  PERCENT 

-  DB GRADED  25  PERCENT 


HORIZONTAL  POSITION 


FIGURE  NO.  13 


a/ 


■  t  •  n 


HORIZONTAL  POSITION 


FIGURE  Nn.  16 
23 


,.iin 


Mil 


200.0 


TIP  DISPL.  VS.  TIME 


1 


w 

Q 

O 

55 


04 

CO 

*-( 

Q 


I 


o 

►— I 

E-* 

« 


■j 


175.0 

150.0 

125.0 

100.0 

75.0 

50.0 

25.0 

0.0 

-25.0 

-50.0 

-75.0 

-100.0 

-125.0 


0.0  20.0  40.0  60.0  80.0  100.0 

TIhE  STEP 


FIGURE  NO.  18 


.  9 

R 


..  w  - 


25 


REFERENCES 


i 


i 

•  j 


i 


I 


vC 


■p 


v, 


* 

ifc 


1.  llerzberg  R.J.,  Johansen  K.F.,  and  Scroud  R.C.,  "Dynamics  and  Control 
of  Large  Satellites,"  Astronautics  and  Aeronautics,  Vol.  16,  Oct., 

1978  pp.  35-39. 

2.  Skelton  R.E.,  "Algorithm  Development  for  the  Control  Design  of  Flexible 
Structures,"  NASA-CP-2258,  1982. 

3.  Skelton,  R.E.,  "Large  Space  System  Control  echnology  Model  Order  Reduction 
Study,"  NASA-CP-2118,  1979. 

4.  Hroner,  C.C.,  "Optimum  Damper  Location  for  a  Free-Free  Beam,"  NASA- 
CP-2168. 

5.  Highsmith,  A.L.,  Stinchcomb,  W.W.,  and  Reifsnider,  K.L.,  "Stiffness 
Reduction  resulting  from  Transverse  cracking  in  Fiber-Reinforced 
Composite  Laminates,"  VPI  &  SU  Report  Mo.  VPI-E-81-33,  Nov.  1981. 

6.  Semi-Annual  Progress  Report  No.  3,  "Cumulative  Damage  Model  for  Advanced 
Composite  Materials,"  General  Dynamics  Technical  Report  FZM-7070, 

October  1982. 

7.  Chou,  P.C.,  Wang,  A.S.D.  and  Miller,  H.,  "Cumulative  Damage  Model  for 
Advanced  Composite  Materials,"  Material  Laboratory  Air  Force  Wright 
Aeronautical  Laboratories,  Report  No.  AFWAL-TR-82-4089 . 

8.  Schapery  R.A. ,  "On  Viscoelastic  deformation  and  failure  behavior  of 
Composite  Materials  with  Distributed  Flaws,"  Advances  in  Aerospace 
Structures  and  Materials-AD-01 ,  1981. 

9.  Allen,  D.H.,  Groves,  S.E.,  and  Schapery,  R.A.,  "A  Damage  Model  For 
Continuous  Fiber  Composites,"  in  preparation. 

10.  Allen,  D.H.,  "A  Prediction  of  Heat  Generation  in  a  Thermoviscoplastic 
Unnxial  Bar,"  Texas  A&M  Mechanics  and  Materials  center,  Report  No.  MM 
4375-83-10. 


11.  Allen,  D.H..  and  Haisler,  W.E . ,  "A  Theory  For  Analysics  of  Thermoplastic 
Materials,"  Computers  and  Structures,  Vol.  13,  1981. 

12.  Bathe,  K.J.,  Wilson,  E.L.,  "Numerical  Methods  in  Finite  Element  Analysis?," 
Prentice-Hall  Inc.,  1976. 

13.  Zienkievicz,  O.C.,  "The  Finite  Element  Method,"  Third  Edition  McGraw- 
Hill  Book  Company,  1978. 

14.  Meiroviteh,  L.,  "  Comnutational  Methods  in  Structural  Dynamics,  Sijthoff 
&  Noordhof,  1480. 


i 


27 


P-«T  W  *\.  ^ 


'W\ 

lbL'Ca?^ 


*v 


,  v' 


PREDICTED  AXIAL  TEMPERATURE  GRADIENT 
IN  A  COUPLED  THERMOVISCOPLASTIC 
UNIAXIAL  BAR 

(WORK  IN  PROGRESS) 


D.H.  Allen 


MM  4875-84-15 


June  1984 


SECURITY  CLASSIFICATION  of  THIS  PAGE  (When  Dml*  Entered) 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


report  number  12.  govt  accession  no. I  3.  Recipient's  catalog  number 


REPORT  DOCUMENTATION  PAGE 


report  number 

MM  4875-84-15 


4.  TITLE  (end  Subtitle) 

Predicted  Axial  Temperature  Gradient  in  a 
Coupled  Thermoviscop las  tic  Uniaxial  Bar 


s.  type  of  report  a  period  covered 


Interim 


.  PERFORMING  ORG.  REPORT  NUMBER 


AUTHOR!"*; 


8.  CONTRACT  OR  GRANT  NUMBER!"*; 


David  H.  Allen 


F49620-83-C-0067 


9.  PERFORMING  ORGANIZATION  NAME  and  AOORESS  ,0-  PROGRAM  element.  PROJECT,  TASK 

Ant  A  a  UNIT  NUMotHs 

Aerospace  Engineering  Department 
Texas  A&M  University 

College  Station,  Texas  77843  _ 


II.  CONTROLLING  OFFICE  NAME  ANO  ADDRESS 

Air  Force  Office  of  Scientific  Research 
Bolling  AFB 

Washington,  D.  C.  20332  _ 

14.  MONITORING  AGENCY  name  4  ACCRF.SSI/1  dillerent  trom  Controlling  Oil  ice)  IS.  SECURITY  CLASS,  (ol  Ihle  report) 

Unc lassif ied 

IS*.  DECLASSIFICATION/  DOWNGRADING 
SCHEDULE 


12.  REPORT  DATE 

June,  1984 


13.  NUMBER  OF  PAGES 


is.  Distribution 


17.  DISTRIBUTION  STATEMENT  -of  the  shslrocl  entered  in  Block  20,  if  dillerent  trom  Report) 


19.  KEY  WORDS  ( Continue  on  reverse  side  7  necessary  end  identity  by  Slock  number) 

thermomechanica.1  coupling 
thermoviscoplastic 
thermomechanical  constitution 
finite  element  method 

heat  generation  _ 


20.  AQSTPACT  'Conf/rm®  :n  reverse  *tJe  If  necessary  and  identify  bv  Slock  number) 

The  thermomechanical  response  of  a  uniaxial  bar  with  thermoviscop Last ic 
constitution  is  predicted  herein  using  the  finite  element  method.  After  a  brief 
review  of  the  governing  field  equations,  variational  principles  are  constructed 
for  the  one  dimensional  conservation  of  momentum  and  energy  equations.  These 
equations  are  coupled  in  that  the  temperature  field  affects  the  displacements 
and  vice  versa. 

Due  to  the  differing  physical  nature  of  the  temperature  and  displacements, 
tirst  order  and  second  order  elements  are  utilized  lot  ihc^e  variables,  rccpec- 


EOITION  IF  NOV  55  IS  OBSOLETE 


SECURITY  CLASSIFICATION  of  This  PAGE  (l»»>*n  Dele  Fnleredi 


SECURITY  CLASSIFICATION  Of  THIS  PAGCfW?l«n  D«l»  gn l trad) 


tively.  The  resulting  semi-discretized  equations  are  then  discretized  in  time 
using  finite  differencing.  This  is  accomplished  by  Euler's  method,  which  is 
utilized  due  to  the  stiff  nature  of  the  constitutive  equations. 

The  model  is  utilized  in  conjunction  with  stress-strain  relations  developed 
by  Bodner  and  Partoin  to  predict  the  axial  temperature  field  in  a  bar  subjected 
to  cyclic  mechanical  end  displacements  and  thermal  flux  boundary  conditions, 
it  is  found  that  spacial  and  time  variation  of  the  temperature  field  is  signifi- 
cnatly  affected  by  the  flux  boundary  conditions. 


PREDICTED  AXIAL  TEMPERATURE  GRADIENT 
IN  A  COUPLED  THERMOVISCOPLASTIC 
UNIAXIAL  BAR 

(Work  in  Progress) 
by 


D.H.  Allen 
Assistant  Professor 
Aerospace  Engineering  Department 
Texas  A&M  University 
College  Station,  TX  77843 


MM  4873-84-15 


June  1984 


ABSTRACT 


The  thermomechanical  response  of  a  uniaxial  bar  with  thermoviscoplastic 
constitution  is  predicted  herein  using  the  finite  element  method.  After  a 
brief  review  of  the  governing  field  equations,  variational  principles  are  con¬ 
structed  for  the  one  dimensional  conservation  of  momentum  and  energy  equations. 

These  equations  are  coupled  in  that  the  temperature  field  affects  the  displace¬ 
ments  and  vice  versa. 

Due  to  the  differing  physical  nature  of  the  temperature  and  displacements, 
first  order  and  second  order  elements  are  utilized  for  these  variables,  respec¬ 
tively.  The  resulting  semi-discretized  equations  are  then  discretized  in  time 
using  finite  differencing.  This  is  accomplished  by  Euler's  method,  which  is 
utilized  due  to  the  stiff  nature  of  the  constitutive  equations. 

The  model  is  utilized  in  conjunction  with  stress-strain  relations  devel¬ 
oped  by  Bodner  and  Partom  ‘o  predict  the  axial  temperature  field  in  a  bar  sub¬ 
jected  to  cyclic  mechanical  end  displacements  and  thermal  flux  boundary  condi¬ 
tions.  It  is  found  that  spacial  and  time  variation  of  the  temperature  field 

i 

is  cignif icantly  affected  by  the  flux  boundary  conditions.  j 

TABLE  OF  SYMBOLS 

time 

I 

axial  internal  resultant  force 
axial  externally  applied  force  per  unit  length 

axial  coordinate  dimension  j 

axial  stress  component 
cross-sectional  area 

end  traction  in  units  of  force  per  unit  area  i 

surface  area 
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S  -  area  of  the  longitudinal  surface  of  the  bar 
£  -  axial  strain  component 

u  -  axial  displacement  component 

ot^  -  internal  state  variable  representing  axial  inelastic  strain 
E  -  Young's  modulus  in  Lie  axial  coordinate  direction 

a  -  coefficient  of  thermal  expansion  in  the  axial  coordinate  direction 
T  -  temperature 

Tg  -  reference  temperature  at  which  no  deformation  is  observed  at  zero  load 

-  internal  state  variable  representing  drag  stress 

q  -  axial  component  heat  flux 

k  -  coefficient  of  axial  thermal  conductivity 

C  -  specific  heat  at  constant  elastic  strain 
v 

P  -  mass  density 

r  -  internal  heat  source  per  unit  mass 
L  -  length  of  the  bar 

INTRODUCTION 

It  is  well  known  that  mechanical  and  thermodynamic  coupling  are  signif¬ 
icant  in  metallic  solids  [1-11].  The  author  has  recently  developed  a  model 
capable  of  predicting  this  coupling  effect  in  thermoviscoplastic  metals  [12]. 
In  the  previous  paper  a  cyclic  strain  control  loading  on  a  sample  of  IN100  at 
1005°K  (1350°F)  was  used  to  predict  a  temperature  rise  of  approximately  3.7°K 
per  cycle  when  the  strain  amplitude  was  2%  and  the  specimen  was  adiabatically 
insulated . 

The  focus  of  the  current  research  is  to  consider  the  effect  of  thermal 
flux  boundary  conditions  on  this  same  process.  The  introduction  of  these  flux 


conditions  causes  the  strain  and  temperature  fields  to  be  inhomogeneous  even 
though  the  stress  field  is  homogeneous  if  the  bar  is  prismatic.  This  spacial 
variation  in  the  field  varables  causes  the  process  to  be  difficult  to  model 
because  the  thermomechanical  constitutive  equations  are  highly  nonlinear  stiff 
differential  equations.  In  this  paper  the  finite  element  method  is  utilized 
to  spatially  discretize  the  dependent  variables,  displacement  and  temperature, 
and  the  finite  difference  method  is  employed  for  timewise  discretization. 

This  process  results  in  a  set  of  highly  nonlinear  algebraic  equations. 

Since  the  thrust  of  this  research  is  to  obtain  accurate  results  without 
regard  to  numerical  efficiency,  the  results  are  obtained  via  the  relatively 
inefficient  but  accurate  method  of  simply  utilizing  successively  smaller  time 
steps  along  with  refined  spatial  mesh  to  obtain  a  convergent  and  therefore 
accurate  solution  for  the  temperature  and  displacement  fields  both  spatially 
and  as  a  function  of  time  for  a  cyclically  imposed  end  displacement. 

The  physical  interest  in  the  problem  is  to  determine  the  effect  of  vari¬ 
ous  flex  boundary  conditions  on  the  predicted  temperature  rise  in  a  bar  sub¬ 
jected  to  cyclic  mechanical  loading.  It  is  found  from  the  analysis  that  the 
introduction  of  flux  boundary  conditions  causes  significant  axial  temperature 
gradients.  Since  flux  conditions  cannot  be  avoided  in  experimental  research, 
it  is  concluded  that  experimental  tests  of  this  type  should  be  viewed  with  cau¬ 
tion  when  their  purpose  is  to  construct  constitutive  relations. 

PROBLEM  SOLUTION 

Field  Problem  Description 
The  following  field  equations  are  given: 


a)  equilibrium  [13], 


where  the  axial  resultant  P  is  defined  by 


P 


odA 


and  (2) 


p  =  i  T  ds  ; 
S 


b)  strain-displacement  relation 


c)  thermomechanical  constitution, 


a=E[e-a1~a(T-  TR) ] 


3a, 

T,  a£)  ,  k  -  l.z 


(3) 


(4) 


(5) 

and  (6) 


q 


(7) 


where  z  is  the  total  number  of  internal  state  variables;  and 

d)  conservation  of  energy 

[<EC  -  ECj  +  Ec,Tr)  +  I  T  £  -  pCy  |I  -  ||  +  pr  -  0  .  (8) 

The  conservation  of  mass  is  satisfied  trivially  and  the  second  law  of  thermo¬ 
dynamics  has  been  previously  shown  to  be  satisfied  by  the  above  equations  [ 14- 
lb  ].  It  should  be  noted  that  equilibrium  equation  (1)  satisfies  equilibrium 
in  the  axial  coordinate  direction  only  in  an  average  sense  over  the  cross-section 


The  above  6+k  equations  (excluding  definition  (3))  define  a  nonlinear 
initial-boundary  value  problem  (together  with  appropriate  thermal  and  mech¬ 
anical  initial  and  boundary  conditions)  in  which  the  following  dependent  vari¬ 
ables  are  sought  as  functions  of  x  and  t:  o,  e,  u,  q,  T,  P,  and  cxR. 

For  convenience  the  domain  is  defined  to  be  of  length  L,  so  that  boundary 
and  initial  conditions  are  of  the  form: 


X 

u(x,0)  =  Uq  =  known 
T  (x ,  0)  ^  T' *  =  known 


litial  conditions  ; 


u(0,t)  =  u^  =  known  or  P(o,t) 


essential  I  u(L,t)  =  u^  =  known  or  P(L,t)  = 
boundary  \ 

conditions  1  T(0,t)  5  T®  *  known  c.  q(0,t)  = 


T(L,t)  = 


Tt  =  known  or  q(L,t) 


P^  =  known 


P  =  known 


qt  =  known 


q^  =  known 


natural 
boundary 
conditions.  (10) 


It  is  now  assumed  that  o  =  a(x)  so  that  equation  (2)  reduces  to 


P  =  OA 


Therefore,  substituting  (4)  into  (5)  and  this  result  into  (11)  gives 


P  =  FA 


r  3u 

[ax  -  ai 


-  oc(T  -  Tr) 


The  above  result  is  now  substituted  into  (1)  to  obtain 


3  f  .  T 3u 

3x  1  LA  Lax  “  ai 


«<T  -  v] } 


-PXU) 


,*•  m'»  '*  •  *»  *  n’V-  *"■  .*•  *’•  »"’•  «,*  -  «.'*  j,*-  .•»**►*<**  *  *.  *  •  *  *  *•  w  * 

*  -  *  «  „  r,  ,  *  .  ■»  «  i  *  «  v*  ^  t  -  *  «  *  «  •  •  ■  •  -_r  *  *f  •  ^  -  -v '  •  ">  mj  -v.  >  *,  i.  \  4 


which  represents  the  differential  equation  relating  displacements  and  temper¬ 


ature  to  the  applied  load  p^(x). 

Equations  (4)  and  (7)  are  next  substituted  into  energy  balance  law  (8) 
and  this  result  is  integrated  over  the  cross-sectional  area  A  to  obtain 


(14) 


where  it  has  been  assumed  that  all  field  variables  depend  on  x  and  t  only. 
Now  define 


(15) 


Careful  inspection  of  equations  (13)  and  (14)  will  indicate  that  these 
equations,  together  with  internal  state  variable  growth  laws  (6)  and  initial 
and  boundary  conditions  (9)  and  (10)  represent  a  well-posed  boundary  value  prob¬ 
lem  in  terms  of  the  2+k  dependent  variables  u,  T,  and  a,. 

Solution  Procedure 

The  field  problem  is  to  be  solved  analytically  using  the  semi-discretized 
rinite  element  technique  with  timewise  finite  differencing.  In  order  to  ac¬ 
complish  this,  differential  equations  (13)  and  (14)  must  first  be  written  in 
a  suitable  variational  form. 

Variational  Equat ions 

Consider  first  equation  (13).  This  governing  equation  is  integrated  against 
a  suitably  smooth  test  function  v  =  v(x)  over  the  domain  of  some  element 


Integrating  by  parts  results  in 


e 


Substituting  equation  (12)  into  the  boundary  term  thus  results  in 


3u 

3x 


a. 


-  a(T  - 


dx  = 


P;VH> 


v(xe) 


P<xe) 


■r 


p  dx 

X 


(18) 


Now  consider  equation  (14).  Once  again  the  governing  equation  is  inte¬ 
grated  against  a  suitably  smooth  test  function  w  =  w(x)  over  the  domain  of  the 

element  £2  : 

e 


(19) 


Integrating  the  heat  flux  term  by  parts  results  in 


-kA  +  wA 

OX  dX 


(E  i  -  Eai  +  e“te)^t-  +  E“21  H 


’EaT  fdx  =  W(xe+1} 


w(xe)  Q(xg) 


/  _  3T 
(PCv  H  - 


pr  )  dx 


where  equation  (15)  has  been  substituted  into  the  boundary  terms. 


Finite  Element  Spacial  Discretization 

Quadratic  displacement  and  linear  temperature  fields  are  now  chosen  within 
each  element: 


C  6 

u(x,t)  =  L  u.  U;.  ,  x  <  x  <  x  ,  and 

,,ii  e  e+1 

1=1 


u(x,t)  -  Z  ll  **  ,  xe  <  x  <  xe+1  , 

i=l 


Q  Q  6  0 

where  u^  =  u^(t)  and  =  T^(t)  are  the  nodal  displacements  and  temperatures, 

6  6  6  6 

respectively,  and  ip^  =  \p_  (x)  and  cfg  =  <J>^(x)  are  quadratic  and  linear  shape 
functions,  respectively  [17]. 

Appropriately,  v  and  w  are  endowed  with  the  properties  of  u  and  T,  respec¬ 
tively,  so  that 
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V  =  ip?  i  =  1,2,3 

w  =  4?  i  =  1 ,2  .  (23) 


Substitution  of  equations  (12)  and  (21)  through  (23)  into  variational  principle 


(18)  results  in 


EA 


dtp? 

l 

dx 


a 


1 


-a{  Z  T?  <f?  -  Td  Jldx  =  -  C(*  P(x  . . )  +  COO  P(*J 


J  =  1 


j  TJ  R 


iv  e+l'  e+1'  i  e  v  e' 


/e+1 

Pxdx  ’ 


i  =  1,2,3 


(24) 


The  above  may  be  written  in  the  following  compact  form 


Z  Ke.  u6  +  Z  S6.  T6  =  Fe  ,  i  =  1,2,3, 

j-1  1J  J  j  =  l  1J  J  1 


where 


EA 


dC  dC 
A  j 

dx  dx 


dx 


e 


d*i  e 

EAa  j —  (J) .  dx 
dx  rj 


e 


i  =  1,2,3,  j  =  1,2,3 


i  1,2,3;  j  =  1,2,3 


(25) 


(26) 


i 

i 


(27) 
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EA 


dx 


F. 

1 


-  (-ex,  +  uT  )dx 

1  K 


-P  (X.) 


ti?  P  dx 
1  x 


i  =  1,2,3. 


(28 


Similarly,  substitution  of  equations  (21)  through  (23)  into  equation  (20)  re 
suits  in 


e 


i  =  1,2. 


(29 


Equations  (29)  may  be  written  in  the  following  form: 


3 

I 


2 

l 


— e  e 

S  .  .  T  + 


due  d/ 
m  m 

dt  dx 


dTe 

m 

dt 


pr 


dx 


■-  / 


e+1  3a. 

A(-EaL  +  EaTR)  •—  +  Q(x.) 


i  =  1,2, 


(30; 


where 


*e+l  di|j?  3a. 

K®  =  J  AE  <()®  dx  i  -  1,2;  j  -  1,2,3;  and  (31] 


*e+l  d<})e  dcf,? 

st,-=-f  kAdJiarLdx  •  1  ■  ‘-2;  j  -  1.2.  (32; 


Finite  Difference  Timewise  Discretization 

Time  dependence  in  equations  (6)  and  (30)  is  handled  via  finite  differ¬ 
encing.  Although  higher  order  approximations  may  be  used,  Euler  forward  dif- 

G  6 

ference  approximations  are  now  entered  for  the  time  rate  of  change  of  a^,  T^, 

and  ue. 
m 


3  a; 


IV 

3t 

(x,t) 

(x, 

,t  + 

At)  • 

dTe 

m 

dt 

(t)  5 

(t 

+ 

At) 

-  Te 
m 

duG 

m 

(t)  3 

[ue  (t 

+ 

At) 

e 

-  u 

e 


(33) 


m  =  1 , 2  ,  and  (34) 


m  =  1,2  . 


(35 


Substitution  of  (33)  through  (35)  into  finite  element  equations  (30)  gives 


£ 


•  e+l 


3a, 


/tlTi  UOA. 

A <p*  [-EttjU)  +  EaTR]  (t) 

e 

+  Q(xi) 


The  above  may  be  written  as  follows: 


I  K..U.+  Z  S..T. 


j=l  13  3  j-1  1J  3 


2  2  2 
+  Z  Z  C.  ..  TS  T®  +  Z  D.  .  Te 
k-1  J-1  ljk  J  k  J-1  1J  J 


2  3  2 

+  Z  Z  E..  .  T®  ue  +  Z  G .  .  T? 
k-1  j-1  lkj  k  3  j-1  1J  3 


e  — o 

+  Z  H .  .  T .  =  FT 


i  =  1,2 


i  =  1,2, 


(37) 


4 


where 


/  A<pl  fr  ♦j  K  dx  * 


i  =  1,2;  j  =  1,2;  k  =  1,2  ,  ( 


Dil  5 


e+l 

-  f  < 


[T*(t)(<J>*)2  +  T®(t)  *®<j>2]  dx, 


i  =  1,2  ,  (39) 


2e+l 

~  J  A ,<|»*  [T®(t)4)*<j>2  +  ^(t)^)2!  dx. 


i  -  1,2  ,  (40) 


E.,  .  5 
ikj 


e+1  g  -e 

-  /  , 


i  =  1,2,;  k  =  1,2;  j  =  1,3 


f  ..e  Ea  f  e.  .  .e  ^1  +  e,  . 

J  A(f,i  Ax  [ul(t)  *1  +  u2(t)  *1 


e  3*2 


w- +  u3(t)  *! 


e  3^1 

iir\ dx>  1 


>L  r 

/  i 


EaTe  e  ^1  e  e^2  e  e  ^3 _ 

At  Ui(t^2^r  +  U2(t)^2^r  +  U3(t)^2^r  dx  »  1 


re+l  pc 

E-  J 


i  =  1,2;  j  =  1,2; 


and  (44) 


/trri 

... 


pc 

e  v 


I  Te(t)4>e  )+  pr  dx 

j=i  J  V 


?e+l 

/  A4>?  [  — Eot  ^ 


(t)  +  EetTR]  j—  (t)  dx 


+  0(x.) 


i  =  1,2 


*.-•  *«:  ■:  '-c-vsv/iv: v:  .v :■  a 
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Equations  (37)  may  be  written  equivalently  as  follows: 


2 

ue  +  Z  S*f.  T?  =  F? 
J  j=1  U  3  x 


£  - 0 

where  K  and  F^  are  as  defined  previously,  and 


2 

+  Z 
k  =  l 


(46) 


and  (47) 


=  _  2 

Sa.  =  S?.  +  Z  C .  . .  T.6  +  D .  .  +  G .  .  +  H .  .  .  (48) 

xj  xj  k=1  ilk  k  xj  xj  ij 


The  above  equations  may  be  a  Ljoined  with  equations  (25)  to  obtain  the  following 
set  of  nonlinear  equation  for  each  element. 


5x5 


where  all  nonlinearity  is  contained  in  [S],  {F  },  and  {F  }. 


Global  Assembly  and  Boundary  Conditions 

Global  assembly  is  accomplished  in  the  standard  way  using  the  Boolean 
matrix  [17).  Interelement  continuity  is  guaranteed  by  setting 


+  1\ 


e+I 


=  0 


.  e 
'2 


+  4* 


e+l 

1 


0 


Boundary  conditions  are  implemented  in  the  standard  way: 
boundary  conditions  are  handled  by  placing  one  on  the  diagonal 


and  (50) 

(51) 

1)  essential 
of  the 


/•  .  *  -  *  •  •  »  '  • 
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appropriate  row  and  zeros  off  diagonal  in  the  stiffness  matrix,  and  the  speci¬ 
fied  value  of  the  essential  variable  on  the  right  hand  side;  and  2)  natural 
bounuary  conditions  are  implemented  directly  to  the  right  hand  side. 


Solution  of  the  Nonlinear  Algebraic  System 

Initial  conditions  are  used  for  the  first  time  step.  The  time  step  At 
is  supplied  for  each  load  increment  and  boundary  conditions  are  incremented 
directly  from  supplied  input  functions. 

The  internal  state  variable  is  handled  in  equations  (23)  and  (45)  by 
using  equations  (35).  is  initialized  as  zero.  The  nonlinear  stiffness 

matrix  [S]  is  initialized  using  nodal  temperatures  and  displacements  from  the 
previous  time  step.  The  displacements  and  temperatures  at  time  t  +  At  are 
then  estimated  directly  and  without  iteration  by  utilizing  equations  (49)  for 
very  small  time  steps. 


CONCLUSION 

The  above  algorithm  has  been  implemented  into  a  fortran  computer  program. 

The  program  is  currently  in  the  final  stage  of  debugging,  with  numerical  results 
expected  in  the  near  future.  It  is  expected  that  the  inclusion  of  thermal  flux 
boundary  conditions  will  cause  a  substantial  decrease  in  the  cyclic  load  in¬ 
duced  temperature  rise  reported  in  reference  12.  These  results  will  be  reported 
as  soon  as  they  become  available. 

Future  research  will  deal  with  the  extension  to  space  structural  applications 
by  including  coordinate  transformation  of  the  truss  elements  and  radiation  bound¬ 
ary  conditions  on  element  longitudinal  surfaces. 
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